TPC/AliTPCtrackerMI.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 //  How to use?  -
  24 //  run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
  25 //  run AliTPCFindTracksMI.C macro - to find tracks
  26 //  tracks are written to AliTPCtracks.root file
  27 //  for comparison also seeds are written to the same file - to special branch
  28 //-------------------------------------------------------
  29
  30
  31 /* $Id$ */
  32
  33
  34
  35 #include "Riostream.h"
  36 #include <TClonesArray.h>
  37 #include <TFile.h>
  38 #include <TObjArray.h>
  39 #include <TTree.h>
  40
  41 #include "AliComplexCluster.h"
  42 #include "AliESD.h"
  43 #include "AliHelix.h"
  44 #include "AliRunLoader.h"
  45 #include "AliTPCClustersRow.h"
  46 #include "AliTPCParam.h"
  47 #include "AliTPCclusterMI.h"
  48 #include "AliTPCpolyTrack.h"
  49 #include "AliTPCreco.h"
  50 #include "AliTPCtrackerMI.h"
  51 #include "TStopwatch.h"
  52
  53 #include "AliTPCReconstructor.h"
  54 #include "AliESDkink.h"
  55 //
  56
  57 ClassImp(AliTPCseed)
  58 ClassImp(AliTPCtrackerMI)
  59
  60
  61 class AliTPCFastMath {
  62 public:
  63   AliTPCFastMath();
  64   static Double_t FastAsin(Double_t x);
  65  private:
  66   static Double_t fgFastAsin[20000];  //lookup table for fast asin computation
  67 };
  68
  69 Double_t AliTPCFastMath::fgFastAsin[20000];
  70 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
  71
  72 AliTPCFastMath::AliTPCFastMath(){
  73   //
  74   // initialized lookup table;
  75   for (Int_t i=0;i<10000;i++){
  76     fgFastAsin[2*i] = TMath::ASin(i/10000.);
  77     fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
  78   }
  79 }
  80
  81 Double_t AliTPCFastMath::FastAsin(Double_t x){
  82   //
  83   // return asin using lookup table
  84   if (x>0){
  85     Int_t index = int(x*10000);
  86     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
  87   }
  88   x*=-1;
  89   Int_t index = int(x*10000);
  90   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
  91 }
  92
  93
  94
  95
  96 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
  97   //
  98   //update track information using current cluster - track->fCurrentCluster
  99
 100
 101   AliTPCclusterMI* c =track->fCurrentCluster;
 102   if (accept>0) track->fCurrentClusterIndex1 |=0x8000;  //sign not accepted clusters
 103
 104   UInt_t i = track->fCurrentClusterIndex1;
 105
 106   Int_t sec=(i&0xff000000)>>24;
 107   //Int_t row = (i&0x00ff0000)>>16;
 108   track->fRow=(i&0x00ff0000)>>16;
 109   track->fSector = sec;
 110   //  Int_t index = i&0xFFFF;
 111   if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
 112   track->SetClusterIndex2(track->fRow, i);
 113   //track->fFirstPoint = row;
 114   //if ( track->fLastPoint<row) track->fLastPoint =row;
 115   //  if (track->fRow<0 || track->fRow>160) {
 116   //  printf("problem\n");
 117   //}
 118   if (track->fFirstPoint>track->fRow)
 119     track->fFirstPoint = track->fRow;
 120   if (track->fLastPoint<track->fRow)
 121     track->fLastPoint  = track->fRow;
 122
 123
 124   track->fClusterPointer[track->fRow] = c;
 125   //
 126
 127   Float_t angle2 = track->GetSnp()*track->GetSnp();
 128   angle2 = TMath::Sqrt(angle2/(1-angle2));
 129   //
 130   //SET NEW Track Point
 131   //
 132   //  if (debug)
 133   {
 134     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->fRow));
 135     //
 136     point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
 137     point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
 138     point.SetErrY(sqrt(track->fErrorY2));
 139     point.SetErrZ(sqrt(track->fErrorZ2));
 140     //
 141     point.SetX(track->GetX());
 142     point.SetY(track->GetY());
 143     point.SetZ(track->GetZ());
 144     point.SetAngleY(angle2);
 145     point.SetAngleZ(track->GetTgl());
 146     if (point.fIsShared){
 147       track->fErrorY2 *= 4;
 148       track->fErrorZ2 *= 4;
 149     }
 150   }
 151
 152   Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
 153   //
 154   track->fErrorY2 *= 1.3;
 155   track->fErrorY2 += 0.01;
 156   track->fErrorZ2 *= 1.3;
 157   track->fErrorZ2 += 0.005;
 158     //}
 159   if (accept>0) return 0;
 160   if (track->GetNumberOfClusters()%20==0){
 161     //    if (track->fHelixIn){
 162     //  TClonesArray & larr = *(track->fHelixIn);
 163     //  Int_t ihelix = larr.GetEntriesFast();
 164     //  new(larr[ihelix]) AliHelix(*track) ;
 165     //}
 166   }
 167   track->fNoCluster =0;
 168   return track->Update(c,chi2,i);
 169 }
 170
 171
 172
 173 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
 174                                       Float_t cory, Float_t corz)
 175 {
 176   //
 177   // decide according desired precision to accept given
 178   // cluster for tracking
 179   Double_t sy2=ErrY2(seed,cluster)*cory;
 180   Double_t sz2=ErrZ2(seed,cluster)*corz;
 181   //sy2=ErrY2(seed,cluster)*cory;
 182   //sz2=ErrZ2(seed,cluster)*cory;
 183
 184   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
 185   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
 186
 187   Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
 188     (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
 189   Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
 190     (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
 191
 192   Double_t rdistance2  = rdistancey2+rdistancez2;
 193   //Int_t  accept =0;
 194
 195   if (rdistance2>16) return 3;
 196
 197
 198   if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
 199     return 2;  //suspisiouce - will be changed
 200
 201   if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
 202     // strict cut on overlaped cluster
 203     return  2;  //suspisiouce - will be changed
 204
 205   if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
 206        && cluster->GetType()<0){
 207     seed->fNFoundable--;
 208     return 2;
 209   }
 210   return 0;
 211 }
 212
 213
 214
 215
 216 //_____________________________________________________________________________
 217 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
 218 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 219 {
 220   //---------------------------------------------------------------------
 221   // The main TPC tracker constructor
 222   //---------------------------------------------------------------------
 223   fInnerSec=new AliTPCSector[fkNIS];
 224   fOuterSec=new AliTPCSector[fkNOS];
 225
 226   Int_t i;
 227   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
 228   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
 229
 230   fN=0;  fSectors=0;
 231
 232   fSeeds=0;
 233   fNtracks = 0;
 234   fParam = par;
 235   Int_t nrowlow = par->GetNRowLow();
 236   Int_t nrowup = par->GetNRowUp();
 237
 238
 239   for (Int_t i=0;i<nrowlow;i++){
 240     fXRow[i]     = par->GetPadRowRadiiLow(i);
 241     fPadLength[i]= par->GetPadPitchLength(0,i);
 242     fYMax[i]     = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
 243   }
 244
 245
 246   for (Int_t i=0;i<nrowup;i++){
 247     fXRow[i+nrowlow]      = par->GetPadRowRadiiUp(i);
 248     fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
 249     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
 250   }
 251   fSeeds=0;
 252   //
 253   fInput    = 0;
 254   fOutput   = 0;
 255   fSeedTree = 0;
 256   fTreeDebug =0;
 257   fNewIO     =0;
 258   fDebug     =0;
 259   fEvent     =0;
 260 }
 261 //________________________________________________________________________
 262 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
 263   AliTracker(t),
 264   fkNIS(t.fkNIS),
 265   fkNOS(t.fkNOS)
 266 {
 267   //------------------------------------
 268   // dummy copy constructor
 269   //------------------------------------------------------------------
 270 }
 271 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
 272   //------------------------------
 273   // dummy
 274   //--------------------------------------------------------------
 275   return *this;
 276 }
 277 //_____________________________________________________________________________
 278 AliTPCtrackerMI::~AliTPCtrackerMI() {
 279   //------------------------------------------------------------------
 280   // TPC tracker destructor
 281   //------------------------------------------------------------------
 282   delete[] fInnerSec;
 283   delete[] fOuterSec;
 284   if (fSeeds) {
 285     fSeeds->Delete();
 286     delete fSeeds;
 287   }
 288 }
 289
 290 void AliTPCtrackerMI::SetIO()
 291 {
 292   //
 293   fNewIO   =  kTRUE;
 294   fInput   =  AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
 295
 296   fOutput  =  AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
 297   if (fOutput){
 298     AliTPCtrack *iotrack= new AliTPCtrack;
 299     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 300     delete iotrack;
 301   }
 302 }
 303
 304
 305 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
 306 {
 307
 308   // set input
 309   fNewIO = kFALSE;
 310   fInput    = 0;
 311   fOutput   = 0;
 312   fSeedTree = 0;
 313   fTreeDebug =0;
 314   fInput = input;
 315   if (input==0){
 316     return;
 317   }
 318   //set output
 319   fOutput = output;
 320   if (output){
 321     AliTPCtrack *iotrack= new AliTPCtrack;
 322     //    iotrack->fHelixIn   = new TClonesArray("AliHelix");
 323     //iotrack->fHelixOut  = new TClonesArray("AliHelix");
 324     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 325     delete iotrack;
 326   }
 327   if (output && (fDebug&2)){
 328     //write the full seed information if specified in debug mode
 329     //
 330     fSeedTree =  new TTree("Seeds","Seeds");
 331     AliTPCseed * vseed = new AliTPCseed;
 332     //
 333     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
 334     arrtr->ExpandCreateFast(160);
 335     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
 336     //
 337     vseed->fPoints = arrtr;
 338     vseed->fEPoints = arre;
 339     //    vseed->fClusterPoints = arrcl;
 340     fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
 341     delete arrtr;
 342     delete arre;
 343     fTreeDebug = new TTree("trackDebug","trackDebug");
 344     TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
 345     fTreeDebug->Branch("debug",&arrd,32000,99);
 346   }
 347
 348
 349   //set ESD event
 350   fEvent  = event;
 351 }
 352
 353 void AliTPCtrackerMI::FillESD(TObjArray* arr)
 354 {
 355   //
 356   //
 357   //fill esds using updated tracks
 358   if (fEvent){
 359     // write tracks to the event
 360     // store index of the track
 361     Int_t nseed=arr->GetEntriesFast();
 362     //FindKinks(arr,fEvent);
 363     for (Int_t i=0; i<nseed; i++) {
 364       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
 365       if (!pt) continue;
 366       pt->UpdatePoints();
 367       pt->PropagateTo(fParam->GetInnerRadiusLow());
 368
 369       if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
 370         AliESDtrack iotrack;
 371         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 372         iotrack.SetTPCPoints(pt->GetPoints());
 373         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 374         //iotrack.SetTPCindex(i);
 375         fEvent->AddTrack(&iotrack);
 376         continue;
 377       }
 378
 379       if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) {
 380         AliESDtrack iotrack;
 381         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 382         iotrack.SetTPCPoints(pt->GetPoints());
 383         //iotrack.SetTPCindex(i);
 384         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 385         fEvent->AddTrack(&iotrack);
 386         continue;
 387       }
 388       //
 389       // short tracks  - maybe decays
 390
 391       if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.70) {
 392         Int_t found,foundable,shared;
 393         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
 394         if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2)){
 395           AliESDtrack iotrack;
 396           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 397           //iotrack.SetTPCindex(i);
 398           iotrack.SetTPCPoints(pt->GetPoints());
 399           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 400           fEvent->AddTrack(&iotrack);
 401           continue;
 402         }
 403       }
 404
 405       if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.8) {
 406         Int_t found,foundable,shared;
 407         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
 408         if (found<20) continue;
 409         if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
 410         //
 411         AliESDtrack iotrack;
 412         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 413         iotrack.SetTPCPoints(pt->GetPoints());
 414         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 415         //iotrack.SetTPCindex(i);
 416         fEvent->AddTrack(&iotrack);
 417         continue;
 418       }
 419       // short tracks  - secondaties
 420       //
 421       if ( (pt->GetNumberOfClusters()>30) ) {
 422         Int_t found,foundable,shared;
 423         pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
 424         if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
 425           AliESDtrack iotrack;
 426           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 427           iotrack.SetTPCPoints(pt->GetPoints());
 428           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 429           //iotrack.SetTPCindex(i);
 430           fEvent->AddTrack(&iotrack);
 431           continue;
 432         }
 433       }
 434
 435       if ( (pt->GetNumberOfClusters()>15)) {
 436         Int_t found,foundable,shared;
 437         pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
 438         if (found<15) continue;
 439         if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
 440         if (float(found)/float(foundable)<0.8) continue;
 441         //
 442         AliESDtrack iotrack;
 443         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 444         iotrack.SetTPCPoints(pt->GetPoints());
 445         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 446         //iotrack.SetTPCindex(i);
 447         fEvent->AddTrack(&iotrack);
 448         continue;
 449       }
 450     }
 451   }
 452   printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
 453 }
 454
 455 void AliTPCtrackerMI::WriteTracks(TTree * tree)
 456 {
 457   //
 458   // write tracks from seed array to selected tree
 459   //
 460   fOutput  = tree;
 461   if (fOutput){
 462     AliTPCtrack *iotrack= new AliTPCtrack;
 463     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 464   }
 465   WriteTracks();
 466 }
 467
 468 void AliTPCtrackerMI::WriteTracks()
 469 {
 470   //
 471   // write tracks to the given output tree -
 472   // output specified with SetIO routine
 473   if (!fSeeds)  return;
 474   if (!fOutput){
 475     SetIO();
 476   }
 477
 478   if (fOutput){
 479     AliTPCtrack *iotrack= 0;
 480     Int_t nseed=fSeeds->GetEntriesFast();
 481     //for (Int_t i=0; i<nseed; i++) {
 482     //  iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
 483     //  if (iotrack) break;
 484     //}
 485     //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 486     TBranch * br = fOutput->GetBranch("tracks");
 487     br->SetAddress(&iotrack);
 488     //
 489     for (Int_t i=0; i<nseed; i++) {
 490       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
 491       if (!pt) continue;
 492       AliTPCtrack * track = new AliTPCtrack(*pt);
 493       iotrack = track;
 494       pt->fLab2 =i;
 495       //      br->SetAddress(&iotrack);
 496       fOutput->Fill();
 497       delete track;
 498       iotrack =0;
 499     }
 500     //fOutput->GetDirectory()->cd();
 501     //fOutput->Write();
 502   }
 503   // delete iotrack;
 504   //
 505   if (fSeedTree){
 506     //write the full seed information if specified in debug mode
 507
 508     AliTPCseed * vseed = new AliTPCseed;
 509     //
 510     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
 511     arrtr->ExpandCreateFast(160);
 512     //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
 513     //arrcl->ExpandCreateFast(160);
 514     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
 515     //
 516     vseed->fPoints = arrtr;
 517     vseed->fEPoints = arre;
 518     //    vseed->fClusterPoints = arrcl;
 519     //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
 520     TBranch * brseed = fSeedTree->GetBranch("seeds");
 521
 522     Int_t nseed=fSeeds->GetEntriesFast();
 523
 524     for (Int_t i=0; i<nseed; i++) {
 525       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
 526       if (!pt) continue;
 527       pt->fPoints = arrtr;
 528       //      pt->fClusterPoints = arrcl;
 529       pt->fEPoints       = arre;
 530       pt->RebuildSeed();
 531       vseed = pt;
 532       brseed->SetAddress(&vseed);
 533       fSeedTree->Fill();
 534       pt->fPoints  = 0;
 535       pt->fEPoints = 0;
 536       //      pt->fClusterPoints = 0;
 537     }
 538     fSeedTree->Write();
 539     if (fTreeDebug) fTreeDebug->Write();
 540   }
 541
 542 }
 543
 544
 545
 546
 547 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
 548   //
 549   //
 550   //seed->SetErrorY2(0.1);
 551   //return 0.1;
 552   //calculate look-up table at the beginning
 553   static Bool_t  ginit = kFALSE;
 554   static Float_t gnoise1,gnoise2,gnoise3;
 555   static Float_t ggg1[10000];
 556   static Float_t ggg2[10000];
 557   static Float_t ggg3[10000];
 558   static Float_t glandau1[10000];
 559   static Float_t glandau2[10000];
 560   static Float_t glandau3[10000];
 561   //
 562   static Float_t gcor01[500];
 563   static Float_t gcor02[500];
 564   static Float_t gcorp[500];
 565   //
 566
 567   //
 568   if (ginit==kFALSE){
 569     for (Int_t i=1;i<500;i++){
 570       Float_t rsigma = float(i)/100.;
 571       gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
 572       gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
 573       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 574     }
 575
 576     //
 577     for (Int_t i=3;i<10000;i++){
 578       //
 579       //
 580       // inner sector
 581       Float_t amp = float(i);
 582       Float_t padlength =0.75;
 583       gnoise1 = 0.0004/padlength;
 584       Float_t nel     = 0.268*amp;
 585       Float_t nprim   = 0.155*amp;
 586       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 587       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 588       if (glandau1[i]>1) glandau1[i]=1;
 589       glandau1[i]*=padlength*padlength/12.;
 590       //
 591       // outer short
 592       padlength =1.;
 593       gnoise2   = 0.0004/padlength;
 594       nel       = 0.3*amp;
 595       nprim     = 0.133*amp;
 596       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 597       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 598       if (glandau2[i]>1) glandau2[i]=1;
 599       glandau2[i]*=padlength*padlength/12.;
 600       //
 601       //
 602       // outer long
 603       padlength =1.5;
 604       gnoise3   = 0.0004/padlength;
 605       nel       = 0.3*amp;
 606       nprim     = 0.133*amp;
 607       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 608       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 609       if (glandau3[i]>1) glandau3[i]=1;
 610       glandau3[i]*=padlength*padlength/12.;
 611       //
 612     }
 613     ginit = kTRUE;
 614   }
 615   //
 616   //
 617   //
 618   Int_t amp = int(TMath::Abs(cl->GetQ()));
 619   if (amp>9999) {
 620     seed->SetErrorY2(1.);
 621     return 1.;
 622   }
 623   Float_t snoise2;
 624   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 625   Int_t ctype = cl->GetType();
 626   Float_t padlength= GetPadPitchLength(seed->fRow);
 627   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 628   angle2 = angle2/(1-angle2);
 629   //
 630   //cluster "quality"
 631   Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
 632   Float_t res;
 633   //
 634   if (fSectors==fInnerSec){
 635     snoise2 = gnoise1;
 636     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 637     if (ctype==0) res *= gcor01[rsigmay];
 638     if ((ctype>0)){
 639       res+=0.002;
 640       res*= gcorp[rsigmay];
 641     }
 642   }
 643   else {
 644     if (padlength<1.1){
 645       snoise2 = gnoise2;
 646       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 647       if (ctype==0) res *= gcor02[rsigmay];
 648       if ((ctype>0)){
 649         res+=0.002;
 650         res*= gcorp[rsigmay];
 651       }
 652     }
 653     else{
 654       snoise2 = gnoise3;
 655       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 656       if (ctype==0) res *= gcor02[rsigmay];
 657       if ((ctype>0)){
 658         res+=0.002;
 659         res*= gcorp[rsigmay];
 660       }
 661     }
 662   }
 663
 664   if (ctype<0){
 665     res+=0.005;
 666     res*=2.4;  // overestimate error 2 times
 667   }
 668   res+= snoise2;
 669
 670   if (res<2*snoise2)
 671     res = 2*snoise2;
 672
 673   seed->SetErrorY2(res);
 674   return res;
 675
 676
 677 }
 678
 679
 680
 681 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 682   //
 683   //
 684   //seed->SetErrorY2(0.1);
 685   //return 0.1;
 686   //calculate look-up table at the beginning
 687   static Bool_t  ginit = kFALSE;
 688   static Float_t gnoise1,gnoise2,gnoise3;
 689   static Float_t ggg1[10000];
 690   static Float_t ggg2[10000];
 691   static Float_t ggg3[10000];
 692   static Float_t glandau1[10000];
 693   static Float_t glandau2[10000];
 694   static Float_t glandau3[10000];
 695   //
 696   static Float_t gcor01[1000];
 697   static Float_t gcor02[1000];
 698   static Float_t gcorp[1000];
 699   //
 700
 701   //
 702   if (ginit==kFALSE){
 703     for (Int_t i=1;i<1000;i++){
 704       Float_t rsigma = float(i)/100.;
 705       gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
 706       gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
 707       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 708     }
 709
 710     //
 711     for (Int_t i=3;i<10000;i++){
 712       //
 713       //
 714       // inner sector
 715       Float_t amp = float(i);
 716       Float_t padlength =0.75;
 717       gnoise1 = 0.0004/padlength;
 718       Float_t nel     = 0.268*amp;
 719       Float_t nprim   = 0.155*amp;
 720       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 721       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 722       if (glandau1[i]>1) glandau1[i]=1;
 723       glandau1[i]*=padlength*padlength/12.;
 724       //
 725       // outer short
 726       padlength =1.;
 727       gnoise2   = 0.0004/padlength;
 728       nel       = 0.3*amp;
 729       nprim     = 0.133*amp;
 730       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 731       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 732       if (glandau2[i]>1) glandau2[i]=1;
 733       glandau2[i]*=padlength*padlength/12.;
 734       //
 735       //
 736       // outer long
 737       padlength =1.5;
 738       gnoise3   = 0.0004/padlength;
 739       nel       = 0.3*amp;
 740       nprim     = 0.133*amp;
 741       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 742       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 743       if (glandau3[i]>1) glandau3[i]=1;
 744       glandau3[i]*=padlength*padlength/12.;
 745       //
 746     }
 747     ginit = kTRUE;
 748   }
 749   //
 750   //
 751   //
 752   Int_t amp = int(TMath::Abs(cl->GetQ()));
 753   if (amp>9999) {
 754     seed->SetErrorY2(1.);
 755     return 1.;
 756   }
 757   Float_t snoise2;
 758   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 759   Int_t ctype = cl->GetType();
 760   Float_t padlength= GetPadPitchLength(seed->fRow);
 761   //
 762   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 763   //  if (angle2<0.6) angle2 = 0.6;
 764   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
 765   //
 766   //cluster "quality"
 767   Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
 768   Float_t res;
 769   //
 770   if (fSectors==fInnerSec){
 771     snoise2 = gnoise1;
 772     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 773     if (ctype==0) res *= gcor01[rsigmaz];
 774     if ((ctype>0)){
 775       res+=0.002;
 776       res*= gcorp[rsigmaz];
 777     }
 778   }
 779   else {
 780     if (padlength<1.1){
 781       snoise2 = gnoise2;
 782       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 783       if (ctype==0) res *= gcor02[rsigmaz];
 784       if ((ctype>0)){
 785         res+=0.002;
 786         res*= gcorp[rsigmaz];
 787       }
 788     }
 789     else{
 790       snoise2 = gnoise3;
 791       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 792       if (ctype==0) res *= gcor02[rsigmaz];
 793       if ((ctype>0)){
 794         res+=0.002;
 795         res*= gcorp[rsigmaz];
 796       }
 797     }
 798   }
 799
 800   if (ctype<0){
 801     res+=0.002;
 802     res*=1.3;
 803   }
 804   if ((ctype<0) &&amp<70){
 805     res+=0.002;
 806     res*=1.3;
 807   }
 808   res += snoise2;
 809   if (res<2*snoise2)
 810      res = 2*snoise2;
 811   if (res>3) res =3;
 812   seed->SetErrorZ2(res);
 813   return res;
 814 }
 815
 816
 817
 818 /*
 819 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 820   //
 821   //
 822   //seed->SetErrorZ2(0.1);
 823   //return 0.1;
 824
 825   Float_t snoise2;
 826   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 827   //
 828   Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
 829   Int_t ctype = cl->GetType();
 830   Float_t amp = TMath::Abs(cl->GetQ());
 831
 832   Float_t nel;
 833   Float_t nprim;
 834   //
 835   Float_t landau=2 ;    //landau fluctuation part
 836   Float_t gg=2;         // gg fluctuation part
 837   Float_t padlength= GetPadPitchLength(seed->GetX());
 838
 839   if (fSectors==fInnerSec){
 840     snoise2 = 0.0004/padlength;
 841     nel     = 0.268*amp;
 842     nprim   = 0.155*amp;
 843     gg      = (2+0.001*nel/(padlength*padlength))/nel;
 844     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 845     if (landau>1) landau=1;
 846   }
 847   else {
 848     snoise2 = 0.0004/padlength;
 849     nel     = 0.3*amp;
 850     nprim   = 0.133*amp;
 851     gg      = (2+0.0008*nel/(padlength*padlength))/nel;
 852     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 853     if (landau>1) landau=1;
 854   }
 855   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
 856
 857   //
 858   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 859   angle2 = TMath::Sqrt((1-angle2));
 860   if (angle2<0.6) angle2 = 0.6;
 861   //angle2 = 1;
 862
 863   Float_t angle = seed->GetTgl()/angle2;
 864   Float_t angular = landau*angle*angle*padlength*padlength/12.;
 865   Float_t res = sdiff + angular;
 866
 867
 868   if ((ctype==0) && (fSectors ==fOuterSec))
 869     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
 870
 871   if ((ctype==0) && (fSectors ==fInnerSec))
 872     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
 873
 874   if ((ctype>0)){
 875     res+=0.005;
 876     res*= TMath::Power(rsigmaz+0.5,1.5);  //0.31+0.147*ctype;
 877   }
 878   if (ctype<0){
 879     res+=0.002;
 880     res*=1.3;
 881   }
 882   if ((ctype<0) &&amp<70){
 883     res+=0.002;
 884     res*=1.3;
 885   }
 886   res += snoise2;
 887   if (res<2*snoise2)
 888      res = 2*snoise2;
 889
 890   seed->SetErrorZ2(res);
 891   return res;
 892 }
 893 */
 894
 895
 896
 897 void AliTPCseed::Reset(Bool_t all)
 898 {
 899   //
 900   //
 901   SetNumberOfClusters(0);
 902   fNFoundable = 0;
 903   SetChi2(0);
 904   ResetCovariance();
 905   /*
 906   if (fTrackPoints){
 907     for (Int_t i=0;i<8;i++){
 908       delete [] fTrackPoints[i];
 909     }
 910     delete fTrackPoints;
 911     fTrackPoints =0;
 912   }
 913   */
 914
 915   if (all){
 916     for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
 917     for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 918   }
 919
 920 }
 921
 922
 923 void AliTPCseed::Modify(Double_t factor)
 924 {
 925
 926   //------------------------------------------------------------------
 927   //This function makes a track forget its history :)
 928   //------------------------------------------------------------------
 929   if (factor<=0) {
 930     ResetCovariance();
 931     return;
 932   }
 933   fC00*=factor;
 934   fC10*=0;  fC11*=factor;
 935   fC20*=0;  fC21*=0;  fC22*=factor;
 936   fC30*=0;  fC31*=0;  fC32*=0;  fC33*=factor;
 937   fC40*=0;  fC41*=0;  fC42*=0;  fC43*=0;  fC44*=factor;
 938   SetNumberOfClusters(0);
 939   fNFoundable =0;
 940   SetChi2(0);
 941   fRemoval = 0;
 942   fCurrentSigmaY2 = 0.000005;
 943   fCurrentSigmaZ2 = 0.000005;
 944   fNoCluster     = 0;
 945   //fFirstPoint = 160;
 946   //fLastPoint  = 0;
 947 }
 948
 949
 950
 951
 952 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
 953 {
 954   //-----------------------------------------------------------------
 955   // This function find proloncation of a track to a reference plane x=xk.
 956   // doesn't change internal state of the track
 957   //-----------------------------------------------------------------
 958
 959   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
 960
 961   if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
 962     return 0;
 963   }
 964
 965   //  Double_t y1=fP0, z1=fP1;
 966   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
 967   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
 968
 969   y = fP0;
 970   z = fP1;
 971   //y += dx*(c1+c2)/(r1+r2);
 972   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 973
 974   Double_t dy = dx*(c1+c2)/(r1+r2);
 975   Double_t dz = 0;
 976   //
 977   Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
 978   /*
 979   if (TMath::Abs(delta)>0.0001){
 980     dz = fP3*TMath::ASin(delta)/fP4;
 981   }else{
 982     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
 983   }
 984   */
 985   dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
 986   //
 987   y+=dy;
 988   z+=dz;
 989
 990
 991   return 1;
 992 }
 993
 994
 995 //_____________________________________________________________________________
 996 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
 997 {
 998   //-----------------------------------------------------------------
 999   // This function calculates a predicted chi2 increment.
1000   //-----------------------------------------------------------------
1001   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
1002   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1003   r00+=fC00; r01+=fC10; r11+=fC11;
1004
1005   Double_t det=r00*r11 - r01*r01;
1006   if (TMath::Abs(det) < 1.e-10) {
1007     Int_t n=GetNumberOfClusters();
1008     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
1009     return 1e10;
1010   }
1011   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1012
1013   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1014
1015   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
1016 }
1017
1018
1019 //_________________________________________________________________________________________
1020
1021
1022 Int_t AliTPCseed::Compare(const TObject *o) const {
1023   //-----------------------------------------------------------------
1024   // This function compares tracks according to the sector - for given sector according z
1025   //-----------------------------------------------------------------
1026   AliTPCseed *t=(AliTPCseed*)o;
1027
1028   if (fSort == 0){
1029     if (t->fRelativeSector>fRelativeSector) return -1;
1030     if (t->fRelativeSector<fRelativeSector) return 1;
1031     Double_t z2 = t->GetZ();
1032     Double_t z1 = GetZ();
1033     if (z2>z1) return 1;
1034     if (z2<z1) return -1;
1035     return 0;
1036   }
1037   else {
1038     Float_t f2 =1;
1039     f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
1040     if (t->fBConstrain) f2=1.2;
1041
1042     Float_t f1 =1;
1043     f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
1044
1045     if (fBConstrain)   f1=1.2;
1046
1047     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
1048     else return +1;
1049   }
1050 }
1051
1052 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
1053 {
1054   //rotate to track "local coordinata
1055   Float_t x = seed->GetX();
1056   Float_t y = seed->GetY();
1057   Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
1058
1059   if (y > ymax) {
1060     seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
1061     if (!seed->Rotate(fSectors->GetAlpha()))
1062       return;
1063   } else if (y <-ymax) {
1064     seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
1065     if (!seed->Rotate(-fSectors->GetAlpha()))
1066       return;
1067   }
1068
1069 }
1070
1071
1072
1073
1074 //_____________________________________________________________________________
1075 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
1076   //-----------------------------------------------------------------
1077   // This function associates a cluster with this track.
1078   //-----------------------------------------------------------------
1079   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1080
1081   r00+=fC00; r01+=fC10; r11+=fC11;
1082   Double_t det=r00*r11 - r01*r01;
1083   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
1084
1085   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
1086   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
1087   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
1088   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
1089   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
1090
1091   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1092   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1093   if (TMath::Abs(cur*fX-eta) >= 0.9) {
1094     return 0;
1095   }
1096
1097   fP0 += k00*dy + k01*dz;
1098   fP1 += k10*dy + k11*dz;
1099   fP2  = eta;
1100   fP3 += k30*dy + k31*dz;
1101   fP4  = cur;
1102
1103   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1104   Double_t c12=fC21, c13=fC31, c14=fC41;
1105
1106   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1107   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
1108   fC40-=k00*c04+k01*c14;
1109
1110   fC11-=k10*c01+k11*fC11;
1111   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
1112   fC41-=k10*c04+k11*c14;
1113
1114   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
1115   fC42-=k20*c04+k21*c14;
1116
1117   fC33-=k30*c03+k31*c13;
1118   fC43-=k40*c03+k41*c13;
1119
1120   fC44-=k40*c04+k41*c14;
1121
1122   Int_t n=GetNumberOfClusters();
1123   //  fIndex[n]=index;
1124   SetNumberOfClusters(n+1);
1125   SetChi2(GetChi2()+chisq);
1126
1127   return 1;
1128 }
1129
1130
1131
1132 //_____________________________________________________________________________
1133 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1134                    Double_t x2,Double_t y2,
1135                    Double_t x3,Double_t y3)
1136 {
1137   //-----------------------------------------------------------------
1138   // Initial approximation of the track curvature
1139   //-----------------------------------------------------------------
1140   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1141   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1142                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1143   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1144                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1145
1146   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1147   if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1148   return -xr*yr/sqrt(xr*xr+yr*yr);
1149 }
1150
1151
1152
1153 //_____________________________________________________________________________
1154 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1155                    Double_t x2,Double_t y2,
1156                    Double_t x3,Double_t y3)
1157 {
1158   //-----------------------------------------------------------------
1159   // Initial approximation of the track curvature
1160   //-----------------------------------------------------------------
1161   x3 -=x1;
1162   x2 -=x1;
1163   y3 -=y1;
1164   y2 -=y1;
1165   //
1166   Double_t det = x3*y2-x2*y3;
1167   if (det==0) {
1168     return 100;
1169   }
1170   //
1171   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1172   Double_t x0 = x3*0.5-y3*u;
1173   Double_t y0 = y3*0.5+x3*u;
1174   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1175   if (det<0) c2*=-1;
1176   return c2;
1177 }
1178
1179
1180 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1181                    Double_t x2,Double_t y2,
1182                    Double_t x3,Double_t y3)
1183 {
1184   //-----------------------------------------------------------------
1185   // Initial approximation of the track curvature
1186   //-----------------------------------------------------------------
1187   x3 -=x1;
1188   x2 -=x1;
1189   y3 -=y1;
1190   y2 -=y1;
1191   //
1192   Double_t det = x3*y2-x2*y3;
1193   if (det==0) {
1194     return 100;
1195   }
1196   //
1197   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1198   Double_t x0 = x3*0.5-y3*u;
1199   Double_t y0 = y3*0.5+x3*u;
1200   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1201   if (det<0) c2*=-1;
1202   x0+=x1;
1203   x0*=c2;
1204   return x0;
1205 }
1206
1207
1208
1209 //_____________________________________________________________________________
1210 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1211                    Double_t x2,Double_t y2,
1212                    Double_t x3,Double_t y3)
1213 {
1214   //-----------------------------------------------------------------
1215   // Initial approximation of the track curvature times center of curvature
1216   //-----------------------------------------------------------------
1217   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1218   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1219                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1220   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1221                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1222
1223   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1224
1225   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1226 }
1227
1228 //_____________________________________________________________________________
1229 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1230                    Double_t x2,Double_t y2,
1231                    Double_t z1,Double_t z2)
1232 {
1233   //-----------------------------------------------------------------
1234   // Initial approximation of the tangent of the track dip angle
1235   //-----------------------------------------------------------------
1236   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1237 }
1238
1239
1240 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1241                    Double_t x2,Double_t y2,
1242                    Double_t z1,Double_t z2, Double_t c)
1243 {
1244   //-----------------------------------------------------------------
1245   // Initial approximation of the tangent of the track dip angle
1246   //-----------------------------------------------------------------
1247
1248   //  Double_t angle1;
1249
1250   //angle1    =  (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1251   //
1252   Double_t d  =  TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1253   if (TMath::Abs(d*c*0.5)>1) return 0;
1254   //  Double_t   angle2    =  TMath::ASin(d*c*0.5);
1255   //  Double_t   angle2    =  AliTPCFastMath::FastAsin(d*c*0.5);
1256   Double_t   angle2    = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1257
1258   angle2  = (z1-z2)*c/(angle2*2.);
1259   return angle2;
1260 }
1261
1262 Bool_t   AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1263 {//-----------------------------------------------------------------
1264   // This function find proloncation of a track to a reference plane x=x2.
1265   //-----------------------------------------------------------------
1266
1267   Double_t dx=x2-x1;
1268
1269   if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1270     return kFALSE;
1271   }
1272
1273   Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1274   Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1275   y = x[0];
1276   z = x[1];
1277
1278   Double_t dy = dx*(c1+c2)/(r1+r2);
1279   Double_t dz = 0;
1280   //
1281   Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1282
1283   if (TMath::Abs(delta)>0.01){
1284     dz = x[3]*TMath::ASin(delta)/x[4];
1285   }else{
1286     dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1287   }
1288
1289   //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1290
1291   y+=dy;
1292   z+=dz;
1293
1294   return kTRUE;
1295 }
1296
1297 Int_t  AliTPCtrackerMI::LoadClusters (TTree *tree)
1298 {
1299   //
1300   //
1301   fInput = tree;
1302   return LoadClusters();
1303 }
1304
1305 Int_t  AliTPCtrackerMI::LoadClusters()
1306 {
1307   //
1308   // load clusters to the memory
1309   AliTPCClustersRow *clrow= new AliTPCClustersRow;
1310   clrow->SetClass("AliTPCclusterMI");
1311   clrow->SetArray(0);
1312   clrow->GetArray()->ExpandCreateFast(10000);
1313   //
1314   //  TTree * tree = fClustersArray.GetTree();
1315
1316   TTree * tree = fInput;
1317   TBranch * br = tree->GetBranch("Segment");
1318   br->SetAddress(&clrow);
1319   //
1320   Int_t j=Int_t(tree->GetEntries());
1321   for (Int_t i=0; i<j; i++) {
1322     br->GetEntry(i);
1323     //
1324     Int_t sec,row;
1325     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1326     //
1327     AliTPCRow * tpcrow=0;
1328     Int_t left=0;
1329     if (sec<fkNIS*2){
1330       tpcrow = &(fInnerSec[sec%fkNIS][row]);
1331       left = sec/fkNIS;
1332     }
1333     else{
1334       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1335       left = (sec-fkNIS*2)/fkNOS;
1336     }
1337     if (left ==0){
1338       tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1339       tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1340       for (Int_t i=0;i<tpcrow->fN1;i++)
1341         tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1342     }
1343     if (left ==1){
1344       tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1345       tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1346       for (Int_t i=0;i<tpcrow->fN2;i++)
1347         tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1348     }
1349   }
1350   //
1351   delete clrow;
1352   LoadOuterSectors();
1353   LoadInnerSectors();
1354   return 0;
1355 }
1356
1357
1358 void AliTPCtrackerMI::UnloadClusters()
1359 {
1360   //
1361   // unload clusters from the memory
1362   //
1363   Int_t nrows = fOuterSec->GetNRows();
1364   for (Int_t sec = 0;sec<fkNOS;sec++)
1365     for (Int_t row = 0;row<nrows;row++){
1366       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1367       //      if (tpcrow){
1368       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1369       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1370       //}
1371       tpcrow->ResetClusters();
1372     }
1373   //
1374   nrows = fInnerSec->GetNRows();
1375   for (Int_t sec = 0;sec<fkNIS;sec++)
1376     for (Int_t row = 0;row<nrows;row++){
1377       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1378       //if (tpcrow){
1379       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1380       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1381       //}
1382       tpcrow->ResetClusters();
1383     }
1384
1385   return ;
1386 }
1387
1388
1389 //_____________________________________________________________________________
1390 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1391   //-----------------------------------------------------------------
1392   // This function fills outer TPC sectors with clusters.
1393   //-----------------------------------------------------------------
1394   Int_t nrows = fOuterSec->GetNRows();
1395   UInt_t index=0;
1396   for (Int_t sec = 0;sec<fkNOS;sec++)
1397     for (Int_t row = 0;row<nrows;row++){
1398       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1399       Int_t sec2 = sec+2*fkNIS;
1400       //left
1401       Int_t ncl = tpcrow->fN1;
1402       while (ncl--) {
1403         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1404         index=(((sec2<<8)+row)<<16)+ncl;
1405         tpcrow->InsertCluster(c,index);
1406       }
1407       //right
1408       ncl = tpcrow->fN2;
1409       while (ncl--) {
1410         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1411         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1412         tpcrow->InsertCluster(c,index);
1413       }
1414       //
1415       // write indexes for fast acces
1416       //
1417       for (Int_t i=0;i<510;i++)
1418         tpcrow->fFastCluster[i]=-1;
1419       for (Int_t i=0;i<tpcrow->GetN();i++){
1420         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1421         tpcrow->fFastCluster[zi]=i;  // write index
1422       }
1423       Int_t last = 0;
1424       for (Int_t i=0;i<510;i++){
1425         if (tpcrow->fFastCluster[i]<0)
1426           tpcrow->fFastCluster[i] = last;
1427         else
1428           last = tpcrow->fFastCluster[i];
1429       }
1430     }
1431   fN=fkNOS;
1432   fSectors=fOuterSec;
1433   return 0;
1434 }
1435
1436
1437 //_____________________________________________________________________________
1438 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1439   //-----------------------------------------------------------------
1440   // This function fills inner TPC sectors with clusters.
1441   //-----------------------------------------------------------------
1442   Int_t nrows = fInnerSec->GetNRows();
1443   UInt_t index=0;
1444   for (Int_t sec = 0;sec<fkNIS;sec++)
1445     for (Int_t row = 0;row<nrows;row++){
1446       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1447       //
1448       //left
1449       Int_t ncl = tpcrow->fN1;
1450       while (ncl--) {
1451         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1452         index=(((sec<<8)+row)<<16)+ncl;
1453         tpcrow->InsertCluster(c,index);
1454       }
1455       //right
1456       ncl = tpcrow->fN2;
1457       while (ncl--) {
1458         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1459         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1460         tpcrow->InsertCluster(c,index);
1461       }
1462       //
1463       // write indexes for fast acces
1464       //
1465       for (Int_t i=0;i<510;i++)
1466         tpcrow->fFastCluster[i]=-1;
1467       for (Int_t i=0;i<tpcrow->GetN();i++){
1468         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1469         tpcrow->fFastCluster[zi]=i;  // write index
1470       }
1471       Int_t last = 0;
1472       for (Int_t i=0;i<510;i++){
1473         if (tpcrow->fFastCluster[i]<0)
1474           tpcrow->fFastCluster[i] = last;
1475         else
1476           last = tpcrow->fFastCluster[i];
1477       }
1478
1479     }
1480
1481   fN=fkNIS;
1482   fSectors=fInnerSec;
1483   return 0;
1484 }
1485
1486
1487
1488 //_________________________________________________________________________
1489 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1490   //--------------------------------------------------------------------
1491   //       Return pointer to a given cluster
1492   //--------------------------------------------------------------------
1493   Int_t sec=(index&0xff000000)>>24;
1494   Int_t row=(index&0x00ff0000)>>16;
1495   Int_t ncl=(index&0x00007fff)>>00;
1496
1497   const AliTPCRow * tpcrow=0;
1498   AliTPCclusterMI * clrow =0;
1499
1500   if (sec<fkNIS*2){
1501     tpcrow = &(fInnerSec[sec%fkNIS][row]);
1502     if (tpcrow==0) return 0;
1503
1504     if (sec<fkNIS) {
1505       if (tpcrow->fN1<=ncl) return 0;
1506       clrow = tpcrow->fClusters1;
1507     }
1508     else {
1509       if (tpcrow->fN2<=ncl) return 0;
1510       clrow = tpcrow->fClusters2;
1511     }
1512   }
1513   else {
1514     tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1515     if (tpcrow==0) return 0;
1516
1517     if (sec-2*fkNIS<fkNOS) {
1518       if (tpcrow->fN1<=ncl) return 0;
1519       clrow = tpcrow->fClusters1;
1520     }
1521     else {
1522       if (tpcrow->fN2<=ncl) return 0;
1523       clrow = tpcrow->fClusters2;
1524     }
1525   }
1526
1527   return &(clrow[ncl]);
1528
1529 }
1530
1531
1532
1533 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1534   //-----------------------------------------------------------------
1535   // This function tries to find a track prolongation to next pad row
1536   //-----------------------------------------------------------------
1537   //
1538   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1539   AliTPCclusterMI *cl=0;
1540   Int_t tpcindex= t.GetClusterIndex2(nr);
1541   //
1542   // update current shape info every 5 pad-row
1543   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1544     GetShape(&t,nr);
1545     //}
1546   //
1547   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters
1548     //
1549     if (tpcindex==-1) return 0; //track in dead zone
1550     if (tpcindex>0){     //
1551       cl = t.fClusterPointer[nr];
1552       if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1553       t.fCurrentClusterIndex1 = tpcindex;
1554     }
1555     if (cl){
1556       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
1557       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1558       //
1559       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1560       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1561
1562       if (TMath::Abs(angle-t.GetAlpha())>0.001){
1563         Double_t rotation = angle-t.GetAlpha();
1564         t.fRelativeSector= relativesector;
1565         t.Rotate(rotation);
1566       }
1567       t.PropagateTo(x);
1568       //
1569       t.fCurrentCluster = cl;
1570       t.fRow = nr;
1571       Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1572       if ((tpcindex&0x8000)==0) accept =0;
1573       if (accept<3) {
1574         //if founded cluster is acceptible
1575         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
1576           t.fErrorY2 += 0.03;
1577           t.fErrorZ2 += 0.03;
1578           t.fErrorY2 *= 3;
1579           t.fErrorZ2 *= 3;
1580         }
1581         t.fNFoundable++;
1582         UpdateTrack(&t,accept);
1583         return 1;
1584       }
1585     }
1586   }
1587   if (fIteration>1) return 0;  // not look for new cluster during refitting
1588   //
1589   UInt_t index=0;
1590   if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;
1591   Double_t  y=t.GetYat(x);
1592   if (TMath::Abs(y)>ymax){
1593     if (y > ymax) {
1594       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1595       if (!t.Rotate(fSectors->GetAlpha()))
1596         return 0;
1597     } else if (y <-ymax) {
1598       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1599       if (!t.Rotate(-fSectors->GetAlpha()))
1600         return 0;
1601     }
1602     //return 1;
1603   }
1604   //
1605   if (!t.PropagateTo(x)) {
1606     if (fIteration==0) t.fRemoval = 10;
1607     return 0;
1608   }
1609   y=t.GetY();
1610   Double_t z=t.GetZ();
1611   //
1612   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1613   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1614   Double_t  roady  =1.;
1615   Double_t  roadz = 1.;
1616   //
1617   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1618     t.fInDead = kTRUE;
1619     t.SetClusterIndex2(nr,-1);
1620     return 0;
1621   }
1622   else
1623     {
1624       if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10)) t.fNFoundable++;
1625       else
1626         return 0;
1627     }
1628   //calculate
1629   if (krow) {
1630     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1631     cl = krow.FindNearest2(y,z,roady,roadz,index);
1632     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1633   }
1634   if (cl) {
1635     t.fCurrentCluster = cl;
1636     t.fRow = nr;
1637     if (fIteration==2&&cl->IsUsed(10)) return 0;
1638     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1639     if (fIteration==2&&cl->IsUsed(11)) {
1640       t.fErrorY2 += 0.03;
1641       t.fErrorZ2 += 0.03;
1642       t.fErrorY2 *= 3;
1643       t.fErrorZ2 *= 3;
1644     }
1645     /*
1646     if (t.fCurrentCluster->IsUsed(10)){
1647       //
1648       //
1649
1650       t.fNShared++;
1651       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1652         t.fRemoval =10;
1653         return 0;
1654       }
1655     }
1656     */
1657     if (accept<3) UpdateTrack(&t,accept);
1658
1659   } else {
1660     if ( fIteration==0 && t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1661
1662   }
1663   return 1;
1664 }
1665
1666 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1667   //-----------------------------------------------------------------
1668   // This function tries to find a track prolongation to next pad row
1669   //-----------------------------------------------------------------
1670   //
1671   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1672   Double_t y,z;
1673   if (!t.GetProlongation(x,y,z)) {
1674     t.fRemoval = 10;
1675     return 0;
1676   }
1677   //
1678   //
1679   if (TMath::Abs(y)>ymax){
1680
1681     if (y > ymax) {
1682       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1683       if (!t.Rotate(fSectors->GetAlpha()))
1684         return 0;
1685     } else if (y <-ymax) {
1686       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1687       if (!t.Rotate(-fSectors->GetAlpha()))
1688         return 0;
1689     }
1690     if (!t.PropagateTo(x)) {
1691       return 0;
1692     }
1693     t.GetProlongation(x,y,z);
1694   }
1695   //
1696   // update current shape info every 3 pad-row
1697   if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1698     //    t.fCurrentSigmaY = GetSigmaY(&t);
1699     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1700     GetShape(&t,nr);
1701   }
1702   //
1703   AliTPCclusterMI *cl=0;
1704   UInt_t index=0;
1705
1706
1707   //Int_t nr2 = nr;
1708   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1709   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1710   Double_t  roady  =1.;
1711   Double_t  roadz = 1.;
1712   //
1713   Int_t row = nr;
1714   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1715     t.fInDead = kTRUE;
1716     t.SetClusterIndex2(row,-1);
1717     return 0;
1718   }
1719   else
1720     {
1721       if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1722     }
1723   //calculate
1724
1725   if ((cl==0)&&(krow)) {
1726     //    cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1727     cl = krow.FindNearest2(y,z,roady,roadz,index);
1728
1729     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1730   }
1731
1732   if (cl) {
1733     t.fCurrentCluster = cl;
1734     //    Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1735     //if (accept<3){
1736       t.SetClusterIndex2(row,index);
1737       t.fClusterPointer[row] = cl;
1738       //}
1739   }
1740   return 1;
1741 }
1742
1743
1744
1745 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
1746   //-----------------------------------------------------------------
1747   // This function tries to find a track prolongation to next pad row
1748   //-----------------------------------------------------------------
1749   t.fCurrentCluster  = 0;
1750   t.fCurrentClusterIndex1 = 0;
1751
1752   Double_t xt=t.GetX();
1753   Int_t     row = GetRowNumber(xt)-1;
1754   Double_t  ymax= GetMaxY(nr);
1755
1756   if (row < nr) return 1; // don't prolongate if not information until now -
1757   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1758     t.fRemoval =10;
1759     return 0;  // not prolongate strongly inclined tracks
1760   }
1761   if (TMath::Abs(t.GetSnp())>0.95) {
1762     t.fRemoval =10;
1763     return 0;  // not prolongate strongly inclined tracks
1764   }
1765
1766   Double_t x= GetXrow(nr);
1767   Double_t y,z;
1768   //t.PropagateTo(x+0.02);
1769   //t.PropagateTo(x+0.01);
1770   if (!t.PropagateTo(x)){
1771     return 0;
1772   }
1773   //
1774   y=t.GetY();
1775   z=t.GetZ();
1776
1777   if (TMath::Abs(y)>ymax){
1778     if (y > ymax) {
1779       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1780       if (!t.Rotate(fSectors->GetAlpha()))
1781         return 0;
1782     } else if (y <-ymax) {
1783       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1784       if (!t.Rotate(-fSectors->GetAlpha()))
1785         return 0;
1786     }
1787     //    if (!t.PropagateTo(x)){
1788     //  return 0;
1789     //}
1790     return 1;
1791     //y = t.GetY();
1792   }
1793   //
1794
1795   AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1796
1797   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1798     t.fInDead = kTRUE;
1799     t.SetClusterIndex2(nr,-1);
1800     return 0;
1801   }
1802   else
1803     {
1804       if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10)) t.fNFoundable++;
1805       else
1806         return 0;
1807     }
1808
1809   // update current
1810   if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1811     //    t.fCurrentSigmaY = GetSigmaY(&t);
1812     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1813     GetShape(&t,nr);
1814   }
1815
1816   AliTPCclusterMI *cl=0;
1817   UInt_t index=0;
1818   //
1819   Double_t roady = 1.;
1820   Double_t roadz = 1.;
1821   //
1822
1823   if (!cl){
1824     index = t.GetClusterIndex2(nr);
1825     if ( (index>0) && (index&0x8000)==0){
1826       cl = t.fClusterPointer[nr];
1827       if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1828       t.fCurrentClusterIndex1 = index;
1829       if (cl) {
1830         t.fCurrentCluster  = cl;
1831         return 1;
1832       }
1833     }
1834   }
1835
1836   if (krow) {
1837     //cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1838     cl = krow.FindNearest2(y,z,roady,roadz,index);
1839   }
1840
1841   if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1842   t.fCurrentCluster  = cl;
1843
1844   return 1;
1845 }
1846
1847
1848 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1849   //-----------------------------------------------------------------
1850   // This function tries to find a track prolongation to next pad row
1851   //-----------------------------------------------------------------
1852
1853   //update error according neighborhoud
1854
1855   if (t.fCurrentCluster) {
1856     t.fRow = nr;
1857     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1858
1859     if (t.fCurrentCluster->IsUsed(10)){
1860       //
1861       //
1862       //  t.fErrorZ2*=2;
1863       //  t.fErrorY2*=2;
1864       t.fNShared++;
1865       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1866         t.fRemoval =10;
1867         return 0;
1868       }
1869     }
1870     if (fIteration>0) accept = 0;
1871     if (accept<3)  UpdateTrack(&t,accept);
1872
1873   } else {
1874     if (fIteration==0){
1875       if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1876       if (  t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1877
1878       if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1879     }
1880   }
1881   return 1;
1882 }
1883
1884
1885
1886 //_____________________________________________________________________________
1887 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1888   //-----------------------------------------------------------------
1889   // This function tries to find a track prolongation.
1890   //-----------------------------------------------------------------
1891   Double_t xt=t.GetX();
1892   //
1893   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1894   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1895   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1896   //
1897   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1898
1899   Int_t first = GetRowNumber(xt)-1;
1900   for (Int_t nr= first; nr>=rf; nr-=step) {
1901     // update kink info
1902     if (t.GetKinkIndexes()[0]>0){
1903       for (Int_t i=0;i<3;i++){
1904         Int_t index = t.GetKinkIndexes()[i];
1905         if (index==0) break;
1906         if (index<0) continue;
1907         //
1908         AliESDkink * kink = fEvent->GetKink(index-1);
1909         if (!kink){
1910           printf("PROBLEM\n");
1911         }
1912         else{
1913           Int_t kinkrow = kink->fRow0+Int_t(1./(0.1+kink->fAngle[2]));
1914           if (kinkrow==nr){
1915             AliExternalTrackParam paramd(t);
1916             kink->SetDaughter(paramd);
1917             kink->Update();
1918             kink->fStatus+=100;
1919           }
1920         }
1921       }
1922     }
1923
1924     if (nr==80) t.UpdateReference();
1925     if (nr<fInnerSec->GetNRows())
1926       fSectors = fInnerSec;
1927     else
1928       fSectors = fOuterSec;
1929     if (FollowToNext(t,nr)==0)
1930       if (!t.IsActive())
1931         return 0;
1932
1933   }
1934   return 1;
1935 }
1936
1937
1938 //_____________________________________________________________________________
1939 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1940   //-----------------------------------------------------------------
1941   // This function tries to find a track prolongation.
1942   //-----------------------------------------------------------------
1943   Double_t xt=t.GetX();
1944   //
1945   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1946   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1947   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1948   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1949
1950   for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1951
1952     if (FollowToNextFast(t,nr)==0)
1953       if (!t.IsActive()) return 0;
1954
1955   }
1956   return 1;
1957 }
1958
1959
1960
1961
1962
1963 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1964   //-----------------------------------------------------------------
1965   // This function tries to find a track prolongation.
1966   //-----------------------------------------------------------------
1967   //  Double_t xt=t.GetX();
1968   //
1969   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1970   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1971   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1972   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1973
1974   Int_t first = t.fFirstPoint;
1975   //
1976   if (first<0) first=0;
1977   for (Int_t nr=first; nr<=rf; nr++) {
1978     //if ( (t.GetSnp()<0.9))
1979     if (t.GetKinkIndexes()[0]<0){
1980       for (Int_t i=0;i<3;i++){
1981         Int_t index = t.GetKinkIndexes()[i];
1982         if (index==0) break;
1983         if (index>0) continue;
1984         index = TMath::Abs(index);
1985         AliESDkink * kink = fEvent->GetKink(index-1);
1986         if (!kink){
1987           printf("PROBLEM\n");
1988         }
1989         else{
1990           Int_t kinkrow = kink->fRow0-Int_t(1./(0.1+kink->fAngle[2]));
1991           if (kinkrow==nr){
1992             AliExternalTrackParam paramm(t);
1993             kink->SetMother(paramm);
1994             kink->Update();
1995             kink->fStatus+=10;
1996           }
1997         }
1998       }
1999     }
2000     if (nr<fInnerSec->GetNRows())
2001       fSectors = fInnerSec;
2002     else
2003       fSectors = fOuterSec;
2004     FollowToNext(t,nr);
2005   }
2006   return 1;
2007 }
2008
2009
2010
2011
2012
2013 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2014 {
2015   //
2016   //
2017   sum1=0;
2018   sum2=0;
2019   Int_t sum=0;
2020   //
2021   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2022   dz2*=dz2;
2023
2024   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2025   dy2*=dy2;
2026   Float_t distance = TMath::Sqrt(dz2+dy2);
2027   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2028
2029   //  Int_t offset =0;
2030   Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
2031   Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
2032   if (lastpoint>160)
2033     lastpoint =160;
2034   if (firstpoint<0)
2035     firstpoint = 0;
2036   if (firstpoint>lastpoint) {
2037     firstpoint =lastpoint;
2038     //    lastpoint  =160;
2039   }
2040
2041
2042   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2043     if (s1->GetClusterIndex2(i)>0) sum1++;
2044     if (s2->GetClusterIndex2(i)>0) sum2++;
2045     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2046       sum++;
2047     }
2048   }
2049   if (sum<5) return 0;
2050
2051   Float_t summin = TMath::Min(sum1+1,sum2+1);
2052   Float_t ratio = (sum+1)/Float_t(summin);
2053   return ratio;
2054 }
2055
2056 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2057 {
2058   //
2059   //
2060   if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
2061   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
2062
2063   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2064   dz2*=dz2;
2065   Float_t dy2 =(s1->GetY() - s2->GetY());
2066   dy2*=dy2;
2067   Float_t distance = dz2+dy2;
2068   if (distance>325.) return ; // if there are far away  - not overlap - to reduce combinatorics
2069
2070   //
2071   Int_t sumshared=0;
2072   //
2073   Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
2074   Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
2075   //
2076   if (firstpoint>=lastpoint-5) return;;
2077
2078   for (Int_t i=firstpoint;i<lastpoint;i++){
2079     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2080     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2081       sumshared++;
2082     }
2083   }
2084   if (sumshared>4){
2085     // sign clusters
2086     //
2087     for (Int_t i=firstpoint;i<lastpoint;i++){
2088       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2089       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2090         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2091         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);;
2092         if (s1->IsActive()&&s2->IsActive()){
2093           p1->fIsShared = kTRUE;
2094           p2->fIsShared = kTRUE;
2095         }
2096       }
2097     }
2098   }
2099   //
2100   if (sumshared>10){
2101     for (Int_t i=0;i<4;i++){
2102       if (s1->fOverlapLabels[3*i]==0){
2103         s1->fOverlapLabels[3*i] = s2->GetLabel();
2104         s1->fOverlapLabels[3*i+1] = sumshared;
2105         s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
2106         break;
2107       }
2108     }
2109     for (Int_t i=0;i<4;i++){
2110       if (s2->fOverlapLabels[3*i]==0){
2111         s2->fOverlapLabels[3*i] = s1->GetLabel();
2112         s2->fOverlapLabels[3*i+1] = sumshared;
2113         s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
2114         break;
2115       }
2116     }
2117   }
2118
2119 }
2120
2121 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
2122 {
2123   //
2124   //sort trackss according sectors
2125   //
2126   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2127     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2128     if (!pt) continue;
2129     //if (pt) RotateToLocal(pt);
2130     pt->fSort = 0;
2131   }
2132   arr->UnSort();
2133   arr->Sort();  // sorting according z
2134   arr->Expand(arr->GetEntries());
2135   //
2136   //
2137   Int_t nseed=arr->GetEntriesFast();
2138   for (Int_t i=0; i<nseed; i++) {
2139     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2140     if (!pt) continue;
2141     for (Int_t j=0;j<=12;j++){
2142       pt->fOverlapLabels[j] =0;
2143     }
2144   }
2145   for (Int_t i=0; i<nseed; i++) {
2146     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2147     if (!pt) continue;
2148     if (pt->fRemoval>10) continue;
2149     for (Int_t j=i+1; j<nseed; j++){
2150       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2151       //      if (pt2){
2152       if (pt2->fRemoval<=10) {
2153         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2154         SignShared(pt,pt2);
2155       }
2156     }
2157   }
2158 }
2159
2160 void  AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2,  Int_t removalindex)
2161 {
2162   //
2163   //sort trackss according sectors
2164   //
2165   if (fDebug&1) {
2166     Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries());
2167   }
2168   //
2169   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2170     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2171     if (!pt) continue;
2172     pt->fSort = 0;
2173   }
2174   arr->UnSort();
2175   arr->Sort();  // sorting according z
2176   arr->Expand(arr->GetEntries());
2177   //
2178   //reset overlap labels
2179   //
2180   Int_t nseed=arr->GetEntriesFast();
2181   for (Int_t i=0; i<nseed; i++) {
2182     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2183     if (!pt) continue;
2184     pt->SetUniqueID(i);
2185     for (Int_t j=0;j<=12;j++){
2186       pt->fOverlapLabels[j] =0;
2187     }
2188   }
2189   //
2190   //sign shared tracks
2191   for (Int_t i=0; i<nseed; i++) {
2192     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2193     if (!pt) continue;
2194     if (pt->fRemoval>10) continue;
2195     Float_t deltac = pt->GetC()*0.1;
2196     for (Int_t j=i+1; j<nseed; j++){
2197       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2198       //      if (pt2){
2199       if (pt2->fRemoval<=10) {
2200         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2201         if (TMath::Abs(pt->GetC()  -pt2->GetC())>deltac) continue;
2202         if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2203         //
2204         SignShared(pt,pt2);
2205       }
2206     }
2207   }
2208   //
2209   // remove highly shared tracks
2210   for (Int_t i=0; i<nseed; i++) {
2211     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2212     if (!pt) continue;
2213     if (pt->fRemoval>10) continue;
2214     //
2215     Int_t sumshared =0;
2216     for (Int_t j=0;j<4;j++){
2217       sumshared = pt->fOverlapLabels[j*3+1];
2218     }
2219     Float_t factor = factor1;
2220     if (pt->fRemoval>0) factor = factor2;
2221     if (sumshared/pt->GetNumberOfClusters()>factor){
2222       for (Int_t j=0;j<4;j++){
2223         if (pt->fOverlapLabels[3*j]==0) continue;
2224         if (pt->fOverlapLabels[3*j+1]<5) continue;
2225         if (pt->fRemoval==removalindex) continue;
2226         AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2227         if (!pt2) continue;
2228         if (pt2->GetSigma2C()<pt->GetSigma2C()){
2229           //      pt->fRemoval = removalindex;
2230           delete arr->RemoveAt(i);
2231           break;
2232         }
2233       }
2234     }
2235   }
2236   arr->Compress();
2237   if (fDebug&1) {
2238     Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries());
2239   }
2240 }
2241
2242
2243
2244
2245
2246
2247 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2248 {
2249   //
2250   //sort tracks in array according mode criteria
2251   Int_t nseed = arr->GetEntriesFast();
2252   for (Int_t i=0; i<nseed; i++) {
2253     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2254     if (!pt) {
2255       continue;
2256     }
2257     pt->fSort = mode;
2258   }
2259   arr->UnSort();
2260   arr->Sort();
2261 }
2262
2263 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t removalindex)
2264 {
2265
2266   //Loop over all tracks and remove "overlaps"
2267   //
2268   //
2269   Int_t nseed = arr->GetEntriesFast();
2270   Int_t good =0;
2271
2272   for (Int_t i=0; i<nseed; i++) {
2273     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2274     if (!pt) {
2275       delete arr->RemoveAt(i);
2276     }
2277     else{
2278       pt->fSort =1;
2279       pt->fBSigned = kFALSE;
2280     }
2281   }
2282   arr->Compress();
2283   nseed = arr->GetEntriesFast();
2284   arr->UnSort();
2285   arr->Sort();
2286   //
2287   //unsign used
2288   UnsignClusters();
2289   //
2290   for (Int_t i=0; i<nseed; i++) {
2291     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2292     if (!pt) {
2293       continue;
2294     }
2295     Int_t found,foundable,shared;
2296     if (pt->IsActive())
2297       pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2298     else
2299       pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2300     //
2301     Double_t factor = factor2;
2302     if (pt->fBConstrain) factor = factor1;
2303
2304     if ((Float_t(shared)/Float_t(found))>factor){
2305       pt->Desactivate(removalindex);
2306       continue;
2307     }
2308
2309     good++;
2310     for (Int_t i=0; i<160; i++) {
2311       Int_t index=pt->GetClusterIndex2(i);
2312       if (index<0 || index&0x8000 ) continue;
2313       AliTPCclusterMI *c= pt->fClusterPointer[i];
2314       if (!c) continue;
2315       //      if (!c->IsUsed(10)) c->Use(10);
2316       //if (pt->IsActive())
2317       c->Use(10);
2318       //else
2319       //        c->Use(5);
2320     }
2321
2322   }
2323   fNtracks = good;
2324   if (fDebug>0){
2325     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2326   }
2327 }
2328
2329
2330 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2331 {
2332
2333   //Loop over all tracks and remove "overlaps"
2334   //
2335   //
2336   UnsignClusters();
2337   //
2338   Int_t nseed = arr->GetEntriesFast();
2339   Float_t * quality = new Float_t[nseed];
2340   Int_t   * indexes = new Int_t[nseed];
2341   Int_t good =0;
2342   //
2343   //
2344   for (Int_t i=0; i<nseed; i++) {
2345     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2346     if (!pt){
2347       quality[i]=-1;
2348       continue;
2349     }
2350     pt->UpdatePoints();    //select first last max dens points
2351     Float_t * points = pt->GetPoints();
2352     if (points[3]<0.8) quality[i] =-1;
2353     //
2354     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2355   }
2356   TMath::Sort(nseed,quality,indexes);
2357   //
2358   //
2359   for (Int_t itrack=0; itrack<nseed; itrack++) {
2360     Int_t trackindex = indexes[itrack];
2361     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2362     if (quality[trackindex]<0){
2363       if (pt) {
2364         delete arr->RemoveAt(trackindex);
2365       }
2366       else{
2367         arr->RemoveAt(trackindex);
2368       }
2369       continue;
2370     }
2371     //
2372     Int_t first = Int_t(pt->GetPoints()[0]);
2373     Int_t last  = Int_t(pt->GetPoints()[2]);
2374     Double_t factor = (pt->fBConstrain) ? factor1: factor2;
2375     //
2376     Int_t found,foundable,shared;
2377     pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
2378     Float_t sharedfactor = Float_t(shared)/Float_t(found);
2379     //
2380     if (Float_t(shared)/Float_t(found)>factor){
2381       if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2382       delete arr->RemoveAt(trackindex);
2383       continue;
2384     }
2385
2386     if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
2387       if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2388       delete arr->RemoveAt(trackindex);
2389       continue;
2390     }
2391
2392     good++;
2393     if (sharedfactor>0.4) continue;
2394     for (Int_t i=first; i<last; i++) {
2395       Int_t index=pt->GetClusterIndex2(i);
2396       // if (index<0 || index&0x8000 ) continue;
2397       if (index<0 || index&0x8000 ) continue;
2398       AliTPCclusterMI *c= pt->fClusterPointer[i];
2399       if (!c) continue;
2400       c->Use(10);
2401     }
2402   }
2403   fNtracks = good;
2404   if (fDebug>0){
2405     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2406   }
2407   delete []quality;
2408   delete []indexes;
2409 }
2410
2411 void AliTPCtrackerMI::UnsignClusters()
2412 {
2413   //
2414   // loop over all clusters and unsign them
2415   //
2416
2417   for (Int_t sec=0;sec<fkNIS;sec++){
2418     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2419       AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2420       for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2421         //      if (cl[icl].IsUsed(10))
2422         cl[icl].Use(-1);
2423       cl = fInnerSec[sec][row].fClusters2;
2424       for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2425         //if (cl[icl].IsUsed(10))
2426           cl[icl].Use(-1);
2427     }
2428   }
2429
2430   for (Int_t sec=0;sec<fkNOS;sec++){
2431     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2432       AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2433       for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2434         //if (cl[icl].IsUsed(10))
2435           cl[icl].Use(-1);
2436       cl = fOuterSec[sec][row].fClusters2;
2437       for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2438         //if (cl[icl].IsUsed(10))
2439         cl[icl].Use(-1);
2440     }
2441   }
2442
2443 }
2444
2445
2446
2447 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2448 {
2449   //
2450   //sign clusters to be "used"
2451   //
2452   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2453   // loop over "primaries"
2454
2455   Float_t sumdens=0;
2456   Float_t sumdens2=0;
2457   Float_t sumn   =0;
2458   Float_t sumn2  =0;
2459   Float_t sumchi =0;
2460   Float_t sumchi2 =0;
2461
2462   Float_t sum    =0;
2463
2464   TStopwatch timer;
2465   timer.Start();
2466
2467   Int_t nseed = arr->GetEntriesFast();
2468   for (Int_t i=0; i<nseed; i++) {
2469     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2470     if (!pt) {
2471       continue;
2472     }
2473     if (!(pt->IsActive())) continue;
2474     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2475     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2476       sumdens += dens;
2477       sumdens2+= dens*dens;
2478       sumn    += pt->GetNumberOfClusters();
2479       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2480       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2481       if (chi2>5) chi2=5;
2482       sumchi  +=chi2;
2483       sumchi2 +=chi2*chi2;
2484       sum++;
2485     }
2486   }
2487
2488   Float_t mdensity = 0.9;
2489   Float_t meann    = 130;
2490   Float_t meanchi  = 1;
2491   Float_t sdensity = 0.1;
2492   Float_t smeann    = 10;
2493   Float_t smeanchi  =0.4;
2494
2495
2496   if (sum>20){
2497     mdensity = sumdens/sum;
2498     meann    = sumn/sum;
2499     meanchi  = sumchi/sum;
2500     //
2501     sdensity = sumdens2/sum-mdensity*mdensity;
2502     sdensity = TMath::Sqrt(sdensity);
2503     //
2504     smeann   = sumn2/sum-meann*meann;
2505     smeann   = TMath::Sqrt(smeann);
2506     //
2507     smeanchi = sumchi2/sum - meanchi*meanchi;
2508     smeanchi = TMath::Sqrt(smeanchi);
2509   }
2510
2511
2512   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2513   //
2514   for (Int_t i=0; i<nseed; i++) {
2515     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2516     if (!pt) {
2517       continue;
2518     }
2519     if (pt->fBSigned) continue;
2520     if (pt->fBConstrain) continue;
2521     //if (!(pt->IsActive())) continue;
2522     /*
2523     Int_t found,foundable,shared;
2524     pt->GetClusterStatistic(0,160,found, foundable,shared);
2525     if (shared/float(found)>0.3) {
2526       if (shared/float(found)>0.9 ){
2527         //delete arr->RemoveAt(i);
2528       }
2529       continue;
2530     }
2531     */
2532     Bool_t isok =kFALSE;
2533     if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2534       isok = kTRUE;
2535     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2536       isok =kTRUE;
2537     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2538       isok =kTRUE;
2539     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2540       isok =kTRUE;
2541
2542     if (isok)
2543       for (Int_t i=0; i<160; i++) {
2544         Int_t index=pt->GetClusterIndex2(i);
2545         if (index<0) continue;
2546         AliTPCclusterMI *c= pt->fClusterPointer[i];
2547         if (!c) continue;
2548         //if (!(c->IsUsed(10))) c->Use();
2549         c->Use(10);
2550       }
2551   }
2552
2553
2554   //
2555   Double_t maxchi  = meanchi+2.*smeanchi;
2556
2557   for (Int_t i=0; i<nseed; i++) {
2558     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2559     if (!pt) {
2560       continue;
2561     }
2562     //if (!(pt->IsActive())) continue;
2563     if (pt->fBSigned) continue;
2564     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2565     if (chi>maxchi) continue;
2566
2567     Float_t bfactor=1;
2568     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2569
2570     //sign only tracks with enoug big density at the beginning
2571
2572     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2573
2574
2575     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2576     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2577
2578     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2579     if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2580       minn=0;
2581
2582     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2583       //Int_t noc=pt->GetNumberOfClusters();
2584       pt->fBSigned = kTRUE;
2585       for (Int_t i=0; i<160; i++) {
2586
2587         Int_t index=pt->GetClusterIndex2(i);
2588         if (index<0) continue;
2589         AliTPCclusterMI *c= pt->fClusterPointer[i];
2590         if (!c) continue;
2591         //      if (!(c->IsUsed(10))) c->Use();
2592         c->Use(10);
2593       }
2594     }
2595   }
2596   //  gLastCheck = nseed;
2597   //  arr->Compress();
2598   if (fDebug>0){
2599     timer.Print();
2600   }
2601 }
2602
2603
2604 void  AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2605 {
2606   // stop not active tracks
2607   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2608   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2609   Int_t nseed = arr->GetEntriesFast();
2610   //
2611   for (Int_t i=0; i<nseed; i++) {
2612     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2613     if (!pt) {
2614       continue;
2615     }
2616     if (!(pt->IsActive())) continue;
2617     StopNotActive(pt,row0,th0, th1,th2);
2618   }
2619 }
2620
2621
2622
2623 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2624  Float_t th2) const
2625 {
2626   // stop not active tracks
2627   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2628   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2629   Int_t sumgood1  = 0;
2630   Int_t sumgood2  = 0;
2631   Int_t foundable = 0;
2632   Int_t maxindex = seed->fLastPoint;  //last foundable row
2633   if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2634     seed->Desactivate(10) ;
2635     return;
2636   }
2637
2638   for (Int_t i=row0; i<maxindex; i++){
2639     Int_t index = seed->GetClusterIndex2(i);
2640     if (index!=-1) foundable++;
2641     //if (!c) continue;
2642     if (foundable<=30) sumgood1++;
2643     if (foundable<=50) {
2644       sumgood2++;
2645     }
2646     else{
2647       break;
2648     }
2649   }
2650   if (foundable>=30.){
2651      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2652   }
2653   if (foundable>=50)
2654     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2655 }
2656
2657
2658 Int_t AliTPCtrackerMI::RefitInward(AliESD *event)
2659 {
2660   //
2661   // back propagation of ESD tracks
2662   //
2663   //return 0;
2664   fEvent = event;
2665   ReadSeeds(event,2);
2666   fIteration=2;
2667   //PrepareForProlongation(fSeeds,1);
2668   PropagateForward2(fSeeds);
2669   Int_t ntracks=0;
2670   Int_t nseed = fSeeds->GetEntriesFast();
2671   for (Int_t i=0;i<nseed;i++){
2672     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2673     if (!seed) continue;
2674     seed->PropagateTo(fParam->GetInnerRadiusLow());
2675     seed->UpdatePoints();
2676     AliESDtrack *esd=event->GetTrack(i);
2677     seed->CookdEdx(0.02,0.6);
2678     CookLabel(seed,0.1); //For comparison only
2679     if (seed->GetNumberOfClusters()>15){
2680       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2681       esd->SetTPCPoints(seed->GetPoints());
2682       ntracks++;
2683     }
2684     else{
2685       //printf("problem\n");
2686     }
2687   }
2688   //FindKinks(fSeeds,event);
2689   Info("RefitInward","Number of refitted tracks %d",ntracks);
2690   fEvent =0;
2691   //WriteTracks();
2692   return 0;
2693 }
2694
2695
2696 Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2697 {
2698   //
2699   // back propagation of ESD tracks
2700   //
2701
2702   fEvent = event;
2703   fIteration = 1;
2704   ReadSeeds(event,0);
2705   PropagateBack(fSeeds);
2706   Int_t nseed = fSeeds->GetEntriesFast();
2707   Int_t ntracks=0;
2708   for (Int_t i=0;i<nseed;i++){
2709     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2710     if (!seed) continue;
2711     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
2712     seed->UpdatePoints();
2713     AliESDtrack *esd=event->GetTrack(i);
2714     seed->CookdEdx(0.02,0.6);
2715     CookLabel(seed,0.1); //For comparison only
2716     if (seed->GetNumberOfClusters()>15){
2717       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2718       esd->SetTPCPoints(seed->GetPoints());
2719       ntracks++;
2720     }
2721   }
2722   //FindKinks(fSeeds,event);
2723   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2724   fEvent =0;
2725   //WriteTracks();
2726   return 0;
2727 }
2728
2729
2730 void AliTPCtrackerMI::DeleteSeeds()
2731 {
2732   //
2733   //delete Seeds
2734   Int_t nseed = fSeeds->GetEntriesFast();
2735   for (Int_t i=0;i<nseed;i++){
2736     AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2737     if (seed) delete fSeeds->RemoveAt(i);
2738   }
2739   delete fSeeds;
2740   fSeeds =0;
2741 }
2742
2743 void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
2744 {
2745   //
2746   //read seeds from the event
2747
2748   Int_t nentr=event->GetNumberOfTracks();
2749   if (fDebug>0){
2750     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2751   }
2752   if (fSeeds)
2753     DeleteSeeds();
2754   if (!fSeeds){
2755     fSeeds = new TObjArray(nentr);
2756   }
2757   UnsignClusters();
2758   //  Int_t ntrk=0;
2759   for (Int_t i=0; i<nentr; i++) {
2760     AliESDtrack *esd=event->GetTrack(i);
2761     ULong_t status=esd->GetStatus();
2762     if (!(status&AliESDtrack::kTPCin)) continue;
2763     AliTPCtrack t(*esd);
2764     AliTPCseed *seed=new AliTPCseed(t/*,t.GetAlpha()*/);seed->ResetClusters();
2765     for (Int_t ikink=0;ikink<3;ikink++) seed->GetKinkIndexes()[ikink] = esd->GetKinkIndex(ikink);
2766     if ((status==AliESDtrack::kTPCin)&&(direction==1)) seed->ResetCovariance();
2767     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
2768     if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2769       fSeeds->AddAt(0,i);
2770       delete seed;
2771       continue;
2772     }
2773     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 )  {
2774       Double_t par0[5],par1[5],x;
2775       esd->GetInnerExternalParameters(x,par0);
2776       esd->GetExternalParameters(x,par1);
2777       Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2778       Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2779       Double_t trdchi2=0;
2780       if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2781       //reset covariance if suspicious
2782       if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2783         seed->ResetCovariance();
2784     }
2785
2786     //
2787     //
2788     // rotate to the local coordinate system
2789
2790     fSectors=fInnerSec; fN=fkNIS;
2791
2792     Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2793     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2794     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2795     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2796     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2797     if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2798     if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2799     alpha-=seed->GetAlpha();
2800     if (!seed->Rotate(alpha)) {
2801       delete seed;
2802       continue;
2803     }
2804     seed->fEsd = esd;
2805     //
2806     //seed->PropagateTo(fSectors->GetX(0));
2807     //
2808     //    Int_t index = esd->GetTPCindex();
2809     //AliTPCseed * seed2= (AliTPCseed*)fSeeds->At(index);
2810     //if (direction==2){
2811     //  AliTPCseed * seed2  = ReSeed(seed,0.,0.5,1.);
2812     //  if (seed2) {
2813     //  delete seed;
2814     //  seed = seed2;
2815     //  }
2816     //}
2817     //
2818     // sign clusters
2819     for (Int_t irow=0;irow<160;irow++){
2820       Int_t index = seed->GetClusterIndex2(irow);
2821       if (index>0){
2822         //
2823         AliTPCclusterMI * cl = GetClusterMI(index);
2824         seed->fClusterPointer[irow] = cl;
2825         if (cl){
2826           if ((index & 0x8000)==0){
2827             cl->Use(10);  // accepted cluster
2828           }else{
2829             cl->Use(6);   // close cluster not accepted
2830           }
2831         }else{
2832            Info("ReadSeeds","Not found cluster");
2833         }
2834       }
2835     }
2836     fSeeds->AddAt(seed,i);
2837   }
2838 }
2839
2840
2841
2842 //_____________________________________________________________________________
2843 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
2844                                  Float_t deltay, Int_t ddsec) {
2845   //-----------------------------------------------------------------
2846   // This function creates track seeds.
2847   // SEEDING WITH VERTEX CONSTRAIN
2848   //-----------------------------------------------------------------
2849   // cuts[0]   - fP4 cut
2850   // cuts[1]   - tan(phi)  cut
2851   // cuts[2]   - zvertex cut
2852   // cuts[3]   - fP3 cut
2853   Int_t nin0  = 0;
2854   Int_t nin1  = 0;
2855   Int_t nin2  = 0;
2856   Int_t nin   = 0;
2857   Int_t nout1 = 0;
2858   Int_t nout2 = 0;
2859
2860   Double_t x[5], c[15];
2861   //  Int_t di = i1-i2;
2862   //
2863   AliTPCseed * seed = new AliTPCseed;
2864   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2865   Double_t cs=cos(alpha), sn=sin(alpha);
2866   //
2867   //  Double_t x1 =fOuterSec->GetX(i1);
2868   //Double_t xx2=fOuterSec->GetX(i2);
2869
2870   Double_t x1 =GetXrow(i1);
2871   Double_t xx2=GetXrow(i2);
2872
2873   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2874
2875   Int_t imiddle = (i2+i1)/2;    //middle pad row index
2876   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2877   const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2878   //
2879   Int_t ns =sec;
2880
2881   const AliTPCRow& kr1=GetRow(ns,i1);
2882   Double_t ymax  = GetMaxY(i1)-kr1.fDeadZone-1.5;
2883   Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2884
2885   //
2886   // change cut on curvature if it can't reach this layer
2887   // maximal curvature set to reach it
2888   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2889   if (dvertexmax*0.5*cuts[0]>0.85){
2890     cuts[0] = 0.85/(dvertexmax*0.5+1.);
2891   }
2892   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
2893
2894   //  Int_t ddsec = 1;
2895   if (deltay>0) ddsec = 0;
2896   // loop over clusters
2897   for (Int_t is=0; is < kr1; is++) {
2898     //
2899     if (kr1[is]->IsUsed(10)) continue;
2900     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2901     //if (TMath::Abs(y1)>ymax) continue;
2902
2903     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2904
2905     // find possible directions
2906     Float_t anglez = (z1-z3)/(x1-x3);
2907     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
2908     //
2909     //
2910     //find   rotation angles relative to line given by vertex and point 1
2911     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2912     Double_t dvertex  = TMath::Sqrt(dvertex2);
2913     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
2914     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);
2915
2916     //
2917     // loop over 2 sectors
2918     Int_t dsec1=-ddsec;
2919     Int_t dsec2= ddsec;
2920     if (y1<0)  dsec2= 0;
2921     if (y1>0)  dsec1= 0;
2922
2923     Double_t dddz1=0;  // direction of delta inclination in z axis
2924     Double_t dddz2=0;
2925     if ( (z1-z3)>0)
2926       dddz1 =1;
2927     else
2928       dddz2 =1;
2929     //
2930     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2931       Int_t sec2 = sec + dsec;
2932       //
2933       //      AliTPCRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2934       //AliTPCRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2935       AliTPCRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
2936       AliTPCRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2937       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2938       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2939
2940       // rotation angles to p1-p3
2941       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2942       Double_t x2,   y2,   z2;
2943       //
2944       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2945
2946       //
2947       Double_t dxx0 =  (xx2-x3)*cs13r;
2948       Double_t dyy0 =  (xx2-x3)*sn13r;
2949       for (Int_t js=index1; js < index2; js++) {
2950         const AliTPCclusterMI *kcl = kr2[js];
2951         if (kcl->IsUsed(10)) continue;
2952         //
2953         //calcutate parameters
2954         //
2955         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
2956         // stright track
2957         if (TMath::Abs(yy0)<0.000001) continue;
2958         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
2959         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2960         Double_t r02 = (0.25+y0*y0)*dvertex2;
2961         //curvature (radius) cut
2962         if (r02<r2min) continue;
2963
2964         nin0++;
2965         //
2966         Double_t c0  = 1/TMath::Sqrt(r02);
2967         if (yy0>0) c0*=-1.;
2968
2969
2970         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
2971         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2972         Double_t dfi0   = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2973         Double_t dfi1   = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2974         //
2975         //
2976         Double_t z0  =  kcl->GetZ();
2977         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
2978         if (TMath::Abs(zzzz2-z0)>0.5) continue;
2979         nin1++;
2980         //
2981         Double_t dip    = (z1-z0)*c0/dfi1;
2982         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2983         //
2984         y2 = kcl->GetY();
2985         if (dsec==0){
2986           x2 = xx2;
2987           z2 = kcl->GetZ();
2988         }
2989         else
2990           {
2991             // rotation
2992             z2 = kcl->GetZ();
2993             x2= xx2*cs-y2*sn*dsec;
2994             y2=+xx2*sn*dsec+y2*cs;
2995           }
2996
2997         x[0] = y1;
2998         x[1] = z1;
2999         x[2] = x0;
3000         x[3] = dip;
3001         x[4] = c0;
3002         //
3003         //
3004         // do we have cluster at the middle ?
3005         Double_t ym,zm;
3006         GetProlongation(x1,xm,x,ym,zm);
3007         UInt_t dummy;
3008         AliTPCclusterMI * cm=0;
3009         if (TMath::Abs(ym)-ymaxm<0){
3010           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3011           if ((!cm) || (cm->IsUsed(10))) {
3012             continue;
3013           }
3014         }
3015         else{
3016           // rotate y1 to system 0
3017           // get state vector in rotated system
3018           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
3019           Double_t xr2  =  x0*cs+yr1*sn*dsec;
3020           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3021           //
3022           GetProlongation(xx2,xm,xr,ym,zm);
3023           if (TMath::Abs(ym)-ymaxm<0){
3024             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3025             if ((!cm) || (cm->IsUsed(10))) {
3026               continue;
3027             }
3028           }
3029         }
3030
3031
3032         Double_t dym = 0;
3033         Double_t dzm = 0;
3034         if (cm){
3035           dym = ym - cm->GetY();
3036           dzm = zm - cm->GetZ();
3037         }
3038         nin2++;
3039
3040
3041         //
3042         //
3043         Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3044         Double_t sy2=kcl->GetSigmaY2()*2.,     sz2=kcl->GetSigmaZ2()*2.;
3045         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3046         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3047         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3048
3049         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3050         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3051         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3052         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3053         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3054         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3055
3056         Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3057         Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3058         Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3059         Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3060
3061         c[0]=sy1;
3062         c[1]=0.;       c[2]=sz1;
3063         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3064         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3065                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3066         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3067         c[13]=f30*sy1*f40+f32*sy2*f42;
3068         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3069
3070         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3071
3072         UInt_t index=kr1.GetIndex(is);
3073         AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
3074
3075         track->fIsSeeding = kTRUE;
3076         track->fSeed1 = i1;
3077         track->fSeed2 = i2;
3078         track->fSeedType=3;
3079
3080
3081         //if (dsec==0) {
3082           FollowProlongation(*track, (i1+i2)/2,1);
3083           Int_t foundable,found,shared;
3084           track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3085           if ((found<0.55*foundable)  || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3086             seed->Reset();
3087             seed->~AliTPCseed();
3088             continue;
3089           }
3090           //}
3091
3092         nin++;
3093         FollowProlongation(*track, i2,1);
3094
3095
3096         //Int_t rc = 1;
3097         track->fBConstrain =1;
3098         //      track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3099         track->fLastPoint = i1;  // first cluster in track position
3100         track->fFirstPoint = track->fLastPoint;
3101
3102         if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3103             track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3104             track->fNShared>0.4*track->GetNumberOfClusters() ) {
3105           seed->Reset();
3106           seed->~AliTPCseed();
3107           continue;
3108         }
3109         nout1++;
3110         // Z VERTEX CONDITION
3111         Double_t zv;
3112         zv = track->GetZ()+track->GetTgl()/track->GetC()*
3113           ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3114         if (TMath::Abs(zv-z3)>cuts[2]) {
3115           FollowProlongation(*track, TMath::Max(i2-20,0));
3116           zv = track->GetZ()+track->GetTgl()/track->GetC()*
3117             ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3118           if (TMath::Abs(zv-z3)>cuts[2]){
3119             FollowProlongation(*track, TMath::Max(i2-40,0));
3120             zv = track->GetZ()+track->GetTgl()/track->GetC()*
3121               ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3122             if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
3123               // make seed without constrain
3124               AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3125               FollowProlongation(*track2, i2,1);
3126               track2->fBConstrain = kFALSE;
3127               track2->fSeedType = 1;
3128               arr->AddLast(track2);
3129               seed->Reset();
3130               seed->~AliTPCseed();
3131               continue;
3132             }
3133             else{
3134               seed->Reset();
3135               seed->~AliTPCseed();
3136               continue;
3137
3138             }
3139           }
3140         }
3141
3142         track->fSeedType =0;
3143         arr->AddLast(track);
3144         seed = new AliTPCseed;
3145         nout2++;
3146         // don't consider other combinations
3147         if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
3148           break;
3149       }
3150     }
3151   }
3152   if (fDebug>3){
3153     Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3154   }
3155   delete seed;
3156 }
3157
3158
3159 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3160                                  Float_t deltay) {
3161
3162
3163
3164   //-----------------------------------------------------------------
3165   // This function creates track seeds.
3166   //-----------------------------------------------------------------
3167   // cuts[0]   - fP4 cut
3168   // cuts[1]   - tan(phi)  cut
3169   // cuts[2]   - zvertex cut
3170   // cuts[3]   - fP3 cut
3171
3172
3173   Int_t nin0  = 0;
3174   Int_t nin1  = 0;
3175   Int_t nin2  = 0;
3176   Int_t nin   = 0;
3177   Int_t nout1 = 0;
3178   Int_t nout2 = 0;
3179   Int_t nout3 =0;
3180   Double_t x[5], c[15];
3181   //
3182   // make temporary seed
3183   AliTPCseed * seed = new AliTPCseed;
3184   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3185   //  Double_t cs=cos(alpha), sn=sin(alpha);
3186   //
3187   //
3188
3189   // first 3 padrows
3190   Double_t x1 = GetXrow(i1-1);
3191   const    AliTPCRow& kr1=GetRow(sec,i1-1);
3192   Double_t y1max  = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
3193   //
3194   Double_t x1p = GetXrow(i1);
3195   const    AliTPCRow& kr1p=GetRow(sec,i1);
3196   //
3197   Double_t x1m = GetXrow(i1-2);
3198   const    AliTPCRow& kr1m=GetRow(sec,i1-2);
3199
3200   //
3201   //last 3 padrow for seeding
3202   AliTPCRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
3203   Double_t    x3   =  GetXrow(i1-7);
3204   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3205   //
3206   AliTPCRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
3207   Double_t    x3p   = GetXrow(i1-6);
3208   //
3209   AliTPCRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
3210   Double_t    x3m   = GetXrow(i1-8);
3211
3212   //
3213   //
3214   // middle padrow
3215   Int_t im = i1-4;                           //middle pad row index
3216   Double_t xm         = GetXrow(im);         // radius of middle pad-row
3217   const AliTPCRow& krm=GetRow(sec,im);   //middle pad -row
3218   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3219   //
3220   //
3221   Double_t deltax  = x1-x3;
3222   Double_t dymax   = deltax*cuts[1];
3223   Double_t dzmax   = deltax*cuts[3];
3224   //
3225   // loop over clusters
3226   for (Int_t is=0; is < kr1; is++) {
3227     //
3228     if (kr1[is]->IsUsed(10)) continue;
3229     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3230     //
3231     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
3232     //
3233     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3234     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3235     //
3236     Double_t y3,   z3;
3237     //
3238     //
3239     UInt_t index;
3240     for (Int_t js=index1; js < index2; js++) {
3241       const AliTPCclusterMI *kcl = kr3[js];
3242       if (kcl->IsUsed(10)) continue;
3243       y3 = kcl->GetY();
3244       // apply angular cuts
3245       if (TMath::Abs(y1-y3)>dymax) continue;
3246       x3 = x3;
3247       z3 = kcl->GetZ();
3248       if (TMath::Abs(z1-z3)>dzmax) continue;
3249       //
3250       Double_t angley = (y1-y3)/(x1-x3);
3251       Double_t anglez = (z1-z3)/(x1-x3);
3252       //
3253       Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3254       Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3255       //
3256       Double_t yyym = angley*(xm-x1)+y1;
3257       Double_t zzzm = anglez*(xm-x1)+z1;
3258
3259       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3260       if (!kcm) continue;
3261       if (kcm->IsUsed(10)) continue;
3262
3263       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3264       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3265       //
3266       //
3267       //
3268       Int_t used  =0;
3269       Int_t found =0;
3270       //
3271       // look around first
3272       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3273                                                       anglez*(x1m-x1)+z1,
3274                                                       erry,errz,index);
3275       //
3276       if (kc1m){
3277         found++;
3278         if (kc1m->IsUsed(10)) used++;
3279       }
3280       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3281                                                       anglez*(x1p-x1)+z1,
3282                                                       erry,errz,index);
3283       //
3284       if (kc1p){
3285         found++;
3286         if (kc1p->IsUsed(10)) used++;
3287       }
3288       if (used>1)  continue;
3289       if (found<1) continue;
3290
3291       //
3292       // look around last
3293       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3294                                                       anglez*(x3m-x3)+z3,
3295                                                       erry,errz,index);
3296       //
3297       if (kc3m){
3298         found++;
3299         if (kc3m->IsUsed(10)) used++;
3300       }
3301       else
3302         continue;
3303       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3304                                                       anglez*(x3p-x3)+z3,
3305                                                       erry,errz,index);
3306       //
3307       if (kc3p){
3308         found++;
3309         if (kc3p->IsUsed(10)) used++;
3310       }
3311       else
3312         continue;
3313       if (used>1)  continue;
3314       if (found<3) continue;
3315       //
3316       Double_t x2,y2,z2;
3317       x2 = xm;
3318       y2 = kcm->GetY();
3319       z2 = kcm->GetZ();
3320       //
3321
3322       x[0]=y1;
3323       x[1]=z1;
3324       x[4]=F1(x1,y1,x2,y2,x3,y3);
3325       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3326       nin0++;
3327       //
3328       x[2]=F2(x1,y1,x2,y2,x3,y3);
3329       nin1++;
3330       //
3331       x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3332       //if (TMath::Abs(x[3]) > cuts[3]) continue;
3333       nin2++;
3334       //
3335       //
3336       Double_t sy1=0.1,  sz1=0.1;
3337       Double_t sy2=0.1,  sz2=0.1;
3338       Double_t sy3=0.1,  sy=0.1, sz=0.1;
3339
3340       Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3341       Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3342       Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3343       Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3344       Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3345       Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3346
3347       Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3348       Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3349       Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3350       Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3351
3352       c[0]=sy1;
3353       c[1]=0.;       c[2]=sz1;
3354       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3355       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3356       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3357       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3358       c[13]=f30*sy1*f40+f32*sy2*f42;
3359       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3360
3361       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3362
3363       UInt_t index=kr1.GetIndex(is);
3364       AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3365
3366       track->fIsSeeding = kTRUE;
3367
3368       nin++;
3369       FollowProlongation(*track, i1-7,1);
3370       if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
3371           track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3372         seed->Reset();
3373         seed->~AliTPCseed();
3374         continue;
3375       }
3376       nout1++;
3377       nout2++;
3378       //Int_t rc = 1;
3379       FollowProlongation(*track, i2,1);
3380       track->fBConstrain =0;
3381       track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3382       track->fFirstPoint = track->fLastPoint;
3383
3384       if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3385           track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
3386           track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3387         seed->Reset();
3388         seed->~AliTPCseed();
3389         continue;
3390       }
3391
3392       {
3393         FollowProlongation(*track, TMath::Max(i2-10,0),1);
3394         AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3395         FollowProlongation(*track2, i2,1);
3396         track2->fBConstrain = kFALSE;
3397         track2->fSeedType = 4;
3398         arr->AddLast(track2);
3399         seed->Reset();
3400         seed->~AliTPCseed();
3401       }
3402
3403
3404       //arr->AddLast(track);
3405       //seed = new AliTPCseed;
3406       nout3++;
3407     }
3408   }
3409
3410   if (fDebug>3){
3411     Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3412   }
3413   delete seed;
3414 }
3415
3416
3417 //_____________________________________________________________________________
3418 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3419                                  Float_t deltay, Bool_t /*bconstrain*/) {
3420   //-----------------------------------------------------------------
3421   // This function creates track seeds - without vertex constraint
3422   //-----------------------------------------------------------------
3423   // cuts[0]   - fP4 cut        - not applied
3424   // cuts[1]   - tan(phi)  cut
3425   // cuts[2]   - zvertex cut    - not applied
3426   // cuts[3]   - fP3 cut
3427   Int_t nin0=0;
3428   Int_t nin1=0;
3429   Int_t nin2=0;
3430   Int_t nin3=0;
3431   //  Int_t nin4=0;
3432   //Int_t nin5=0;
3433
3434
3435
3436   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3437   //  Double_t cs=cos(alpha), sn=sin(alpha);
3438   Int_t row0 = (i1+i2)/2;
3439   Int_t drow = (i1-i2)/2;
3440   const AliTPCRow& kr0=fSectors[sec][row0];
3441   AliTPCRow * kr=0;
3442
3443   AliTPCpolyTrack polytrack;
3444   Int_t nclusters=fSectors[sec][row0];
3445   AliTPCseed * seed = new AliTPCseed;
3446
3447   Int_t sumused=0;
3448   Int_t cused=0;
3449   Int_t cnused=0;
3450   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3451     Int_t nfound =0;
3452     Int_t nfoundable =0;
3453     for (Int_t iter =1; iter<2; iter++){   //iterations
3454       const AliTPCRow& krm=fSectors[sec][row0-iter];
3455       const AliTPCRow& krp=fSectors[sec][row0+iter];
3456       const AliTPCclusterMI * cl= kr0[is];
3457
3458       if (cl->IsUsed(10)) {
3459         cused++;
3460       }
3461       else{
3462         cnused++;
3463       }
3464       Double_t x = kr0.GetX();
3465       // Initialization of the polytrack
3466       nfound =0;
3467       nfoundable =0;
3468       polytrack.Reset();
3469       //
3470       Double_t y0= cl->GetY();
3471       Double_t z0= cl->GetZ();
3472       Float_t erry = 0;
3473       Float_t errz = 0;
3474
3475       Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3476       if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;  // seed only at the edge
3477
3478       erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3479       errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3480       polytrack.AddPoint(x,y0,z0,erry, errz);
3481
3482       sumused=0;
3483       if (cl->IsUsed(10)) sumused++;
3484
3485
3486       Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3487       Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3488       //
3489       x = krm.GetX();
3490       AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3491       if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3492         erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3493         errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3494         if (cl1->IsUsed(10))  sumused++;
3495         polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3496       }
3497       //
3498       x = krp.GetX();
3499       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3500       if (cl2) {
3501         erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3502         errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3503         if (cl2->IsUsed(10)) sumused++;
3504         polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3505       }
3506       //
3507       if (sumused>0) continue;
3508       nin0++;
3509       polytrack.UpdateParameters();
3510       // follow polytrack
3511       roadz = 1.2;
3512       roady = 1.2;
3513       //
3514       Double_t yn,zn;
3515       nfoundable = polytrack.GetN();
3516       nfound     = nfoundable;
3517       //
3518       for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3519         Float_t maxdist = 0.8*(1.+3./(ddrow));
3520         for (Int_t delta = -1;delta<=1;delta+=2){
3521           Int_t row = row0+ddrow*delta;
3522           kr = &(fSectors[sec][row]);
3523           Double_t xn = kr->GetX();
3524           Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3525           polytrack.GetFitPoint(xn,yn,zn);
3526           if (TMath::Abs(yn)>ymax) continue;
3527           nfoundable++;
3528           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3529           if (cln) {
3530             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3531             if (dist<maxdist){
3532               /*
3533               erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3534               errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3535               if (cln->IsUsed(10)) {
3536                 //      printf("used\n");
3537                 sumused++;
3538                 erry*=2;
3539                 errz*=2;
3540               }
3541               */
3542               erry=0.1;
3543               errz=0.1;
3544               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3545               nfound++;
3546             }
3547           }
3548         }
3549         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;
3550         polytrack.UpdateParameters();
3551       }
3552     }
3553     if ( (sumused>3) || (sumused>0.5*nfound))  {
3554       //printf("sumused   %d\n",sumused);
3555       continue;
3556     }
3557     nin1++;
3558     Double_t dy,dz;
3559     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3560     AliTPCpolyTrack track2;
3561
3562     polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3563     if (track2.GetN()<0.5*nfoundable) continue;
3564     nin2++;
3565
3566     if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3567       //
3568       // test seed with and without constrain
3569       for (Int_t constrain=0; constrain<=0;constrain++){
3570         // add polytrack candidate
3571
3572         Double_t x[5], c[15];
3573         Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3574         track2.GetBoundaries(x3,x1);
3575         x2 = (x1+x3)/2.;
3576         track2.GetFitPoint(x1,y1,z1);
3577         track2.GetFitPoint(x2,y2,z2);
3578         track2.GetFitPoint(x3,y3,z3);
3579         //
3580         //is track pointing to the vertex ?
3581         Double_t x0,y0,z0;
3582         x0=0;
3583         polytrack.GetFitPoint(x0,y0,z0);
3584
3585         if (constrain) {
3586           x2 = x3;
3587           y2 = y3;
3588           z2 = z3;
3589
3590           x3 = 0;
3591           y3 = 0;
3592           z3 = 0;
3593         }
3594         x[0]=y1;
3595         x[1]=z1;
3596         x[4]=F1(x1,y1,x2,y2,x3,y3);
3597
3598         //      if (TMath::Abs(x[4]) >= cuts[0]) continue;  //
3599         x[2]=F2(x1,y1,x2,y2,x3,y3);
3600
3601         //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3602         //x[3]=F3(x1,y1,x2,y2,z1,z2);
3603         x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3604         //if (TMath::Abs(x[3]) > cuts[3]) continue;
3605
3606
3607         Double_t sy =0.1, sz =0.1;
3608         Double_t sy1=0.02, sz1=0.02;
3609         Double_t sy2=0.02, sz2=0.02;
3610         Double_t sy3=0.02;
3611
3612         if (constrain){
3613           sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3614         }
3615
3616         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3617         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3618         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3619         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3620         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3621         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3622
3623         Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3624         Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3625         Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3626         Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3627
3628
3629         c[0]=sy1;
3630         c[1]=0.;       c[2]=sz1;
3631         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3632         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3633         c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3634         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3635         c[13]=f30*sy1*f40+f32*sy2*f42;
3636         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3637
3638         //Int_t row1 = fSectors->GetRowNumber(x1);
3639         Int_t row1 = GetRowNumber(x1);
3640
3641         UInt_t index=0;
3642         //kr0.GetIndex(is);
3643         AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3644         track->fIsSeeding = kTRUE;
3645         Int_t rc=FollowProlongation(*track, i2);
3646         if (constrain) track->fBConstrain =1;
3647         else
3648           track->fBConstrain =0;
3649         track->fLastPoint = row1+fInnerSec->GetNRows();  // first cluster in track position
3650         track->fFirstPoint = track->fLastPoint;
3651
3652         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3653             track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3654             track->fNShared>0.4*track->GetNumberOfClusters()) {
3655           //delete track;
3656           seed->Reset();
3657           seed->~AliTPCseed();
3658         }
3659         else {
3660           arr->AddLast(track);
3661           seed = new AliTPCseed;
3662         }
3663         nin3++;
3664       }
3665     }  // if accepted seed
3666   }
3667   if (fDebug>3){
3668     Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3669   }
3670   delete seed;
3671 }
3672
3673
3674 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3675 {
3676   //
3677   //
3678   //reseed using track points
3679   Int_t p0 = int(r0*track->GetNumberOfClusters());     // point 0
3680   Int_t p1 = int(r1*track->GetNumberOfClusters());
3681   Int_t p2 = int(r2*track->GetNumberOfClusters());   // last point
3682   Int_t pp2=0;
3683   Double_t  x0[3],x1[3],x2[3];
3684   x0[0]=-1;
3685   x0[0]=-1;
3686   x0[0]=-1;
3687
3688   // find track position at given ratio of the length
3689   Int_t  sec0, sec1, sec2;
3690   sec0=0;
3691   sec1=0;
3692   sec2=0;
3693   Int_t index=-1;
3694   Int_t clindex;
3695   for (Int_t i=0;i<160;i++){
3696     if (track->fClusterPointer[i]){
3697       index++;
3698       AliTPCTrackerPoint   *trpoint =track->GetTrackPoint(i);
3699       if ( (index<p0) || x0[0]<0 ){
3700         if (trpoint->GetX()>1){
3701           clindex = track->GetClusterIndex2(i);
3702           if (clindex>0){
3703             x0[0] = trpoint->GetX();
3704             x0[1] = trpoint->GetY();
3705             x0[2] = trpoint->GetZ();
3706             sec0  = ((clindex&0xff000000)>>24)%18;
3707           }
3708         }
3709       }
3710
3711       if ( (index<p1) &&(trpoint->GetX()>1)){
3712         clindex = track->GetClusterIndex2(i);
3713         if (clindex>0){
3714           x1[0] = trpoint->GetX();
3715           x1[1] = trpoint->GetY();
3716           x1[2] = trpoint->GetZ();
3717           sec1  = ((clindex&0xff000000)>>24)%18;
3718         }
3719       }
3720       if ( (index<p2) &&(trpoint->GetX()>1)){
3721         clindex = track->GetClusterIndex2(i);
3722         if (clindex>0){
3723           x2[0] = trpoint->GetX();
3724           x2[1] = trpoint->GetY();
3725           x2[2] = trpoint->GetZ();
3726           sec2  = ((clindex&0xff000000)>>24)%18;
3727           pp2 = i;
3728         }
3729       }
3730     }
3731   }
3732
3733   Double_t alpha, cs,sn, xx2,yy2;
3734   //
3735   alpha = (sec1-sec2)*fSectors->GetAlpha();
3736   cs = TMath::Cos(alpha);
3737   sn = TMath::Sin(alpha);
3738   xx2= x1[0]*cs-x1[1]*sn;
3739   yy2= x1[0]*sn+x1[1]*cs;
3740   x1[0] = xx2;
3741   x1[1] = yy2;
3742   //
3743   alpha = (sec0-sec2)*fSectors->GetAlpha();
3744   cs = TMath::Cos(alpha);
3745   sn = TMath::Sin(alpha);
3746   xx2= x0[0]*cs-x0[1]*sn;
3747   yy2= x0[0]*sn+x0[1]*cs;
3748   x0[0] = xx2;
3749   x0[1] = yy2;
3750   //
3751   //
3752   //
3753   Double_t x[5],c[15];
3754   //
3755   x[0]=x2[1];
3756   x[1]=x2[2];
3757   x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3758   //  if (x[4]>1) return 0;
3759   x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3760   x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3761   //if (TMath::Abs(x[3]) > 2.2)  return 0;
3762   //if (TMath::Abs(x[2]) > 1.99) return 0;
3763   //
3764   Double_t sy =0.1,  sz =0.1;
3765   //
3766   Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3767   Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3768   Double_t sy3=0.01+track->GetSigmaY2();
3769   //
3770   Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3771   Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3772   Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3773   Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3774   Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3775   Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3776   //
3777   Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3778   Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3779   Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3780   Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3781
3782
3783   c[0]=sy1;
3784   c[1]=0.;       c[2]=sz1;
3785   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3786   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3787   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3788   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3789   c[13]=f30*sy1*f40+f32*sy2*f42;
3790   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3791
3792   //  Int_t row1 = fSectors->GetRowNumber(x2[0]);
3793   AliTPCseed *seed=new  AliTPCseed(0, x, c, x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3794   //  Double_t y0,z0,y1,z1, y2,z2;
3795   //seed->GetProlongation(x0[0],y0,z0);
3796   // seed->GetProlongation(x1[0],y1,z1);
3797   //seed->GetProlongation(x2[0],y2,z2);
3798   //  seed =0;
3799   seed->fLastPoint  = pp2;
3800   seed->fFirstPoint = pp2;
3801
3802
3803   return seed;
3804 }
3805
3806
3807 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3808 {
3809   //
3810   //
3811   //reseed using founded clusters
3812   //
3813   // Find the number of clusters
3814   Int_t nclusters = 0;
3815   for (Int_t irow=0;irow<160;irow++){
3816     if (track->GetClusterIndex(irow)>0) nclusters++;
3817   }
3818   //
3819   Int_t ipos[3];
3820   ipos[0] = TMath::Max(int(r0*nclusters),0);             // point 0 cluster
3821   ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1);   //
3822   ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1);   // last point
3823   //
3824   //
3825   Double_t  xyz[3][3];
3826   Int_t     row[3],sec[3]={0,0,0};
3827   //
3828   // find track row position at given ratio of the length
3829   Int_t index=-1;
3830   for (Int_t irow=0;irow<160;irow++){
3831     if (track->GetClusterIndex2(irow)<0) continue;
3832     index++;
3833     for (Int_t ipoint=0;ipoint<3;ipoint++){
3834       if (index<=ipos[ipoint]) row[ipoint] = irow;
3835     }
3836   }
3837   //
3838   //Get cluster and sector position
3839   for (Int_t ipoint=0;ipoint<3;ipoint++){
3840     Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3841     AliTPCclusterMI * cl = GetClusterMI(clindex);
3842     if (cl==0) {
3843       //Error("Bug\n");
3844       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
3845       return 0;
3846     }
3847     sec[ipoint]     = ((clindex&0xff000000)>>24)%18;
3848     xyz[ipoint][0]  = GetXrow(row[ipoint]);
3849     xyz[ipoint][1]  = cl->GetY();
3850     xyz[ipoint][2]  = cl->GetZ();
3851   }
3852   //
3853   //
3854   // Calculate seed state vector and covariance matrix
3855
3856   Double_t alpha, cs,sn, xx2,yy2;
3857   //
3858   alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3859   cs = TMath::Cos(alpha);
3860   sn = TMath::Sin(alpha);
3861   xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3862   yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3863   xyz[1][0] = xx2;
3864   xyz[1][1] = yy2;
3865   //
3866   alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3867   cs = TMath::Cos(alpha);
3868   sn = TMath::Sin(alpha);
3869   xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3870   yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3871   xyz[0][0] = xx2;
3872   xyz[0][1] = yy2;
3873   //
3874   //
3875   //
3876   Double_t x[5],c[15];
3877   //
3878   x[0]=xyz[2][1];
3879   x[1]=xyz[2][2];
3880   x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3881   x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3882   x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3883   //
3884   Double_t sy =0.1,  sz =0.1;
3885   //
3886   Double_t sy1=0.2, sz1=0.2;
3887   Double_t sy2=0.2, sz2=0.2;
3888   Double_t sy3=0.2;
3889   //
3890   Double_t f40=(F1(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[4])/sy;
3891   Double_t f42=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[4])/sy;
3892   Double_t f43=(F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[4])/sy;
3893   Double_t f20=(F2(xyz[2][0],xyz[2][1]+sy,xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1])-x[2])/sy;
3894   Double_t f22=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1]+sy,xyz[0][0],xyz[0][1])-x[2])/sy;
3895   Double_t f23=(F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]+sy)-x[2])/sy;
3896   //
3897   Double_t f30=(F3(xyz[2][0],xyz[2][1]+sy,xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2])-x[3])/sy;
3898   Double_t f31=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2]+sz,xyz[0][2])-x[3])/sz;
3899   Double_t f32=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1]+sy,xyz[2][2],xyz[0][2])-x[3])/sy;
3900   Double_t f34=(F3(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2]+sz)-x[3])/sz;
3901
3902
3903   c[0]=sy1;
3904   c[1]=0.;       c[2]=sz1;
3905   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3906   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3907   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3908   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3909   c[13]=f30*sy1*f40+f32*sy2*f42;
3910   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3911
3912   //  Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3913   AliTPCseed *seed=new  AliTPCseed(0, x, c, xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3914   seed->fLastPoint  = row[2];
3915   seed->fFirstPoint = row[2];
3916   return seed;
3917 }
3918
3919 void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
3920 {
3921   //
3922   //  find kinks
3923   //
3924   //
3925   TObjArray *kinks= new TObjArray(10000);
3926   Int_t nentries = array->GetEntriesFast();
3927   AliHelix *helixes      = new AliHelix[nentries];
3928   Int_t    *sign         = new Int_t[nentries];
3929   Int_t    *nclusters    = new Int_t[nentries];
3930   Float_t  *alpha        = new Float_t[nentries];
3931   AliESDkink * kink      = new AliESDkink();
3932   Int_t      * usage     = new Int_t[nentries];
3933   Float_t  *zm             = new Float_t[nentries];
3934   Float_t  *fim            = new Float_t[nentries];
3935
3936   //
3937   //
3938   for (Int_t i=0;i<nentries;i++){
3939     sign[i]=0;
3940     usage[i]=0;
3941     AliTPCseed* track = (AliTPCseed*)array->At(i);
3942     if (!track) continue;
3943     track->UpdatePoints();
3944     if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
3945       nclusters[i]=track->GetNumberOfClusters();
3946       alpha[i] = track->GetAlpha();
3947       new (&helixes[i]) AliHelix(*track);
3948       sign[i] = (track->GetC()>0) ? -1:1;
3949       Double_t x,y,z;
3950       x=160;
3951       if (track->GetProlongation(x,y,z)){
3952         zm[i]  = z;
3953         fim[i] = alpha[i]+TMath::ATan2(y,x);
3954       }
3955       else{
3956         zm[i]  = track->GetZ();
3957         fim[i] = alpha[i];
3958       }
3959     }
3960   }
3961   //
3962   //
3963   TStopwatch timer;
3964   timer.Start();
3965   Int_t ncandidates =0;
3966   Int_t nall =0;
3967   Int_t ntracks=0;
3968   Double_t phase[2][2],radius[2];
3969   //
3970   //
3971   for (Int_t i =0;i<nentries;i++){
3972     if (sign[i]==0) continue;
3973     AliTPCseed * track0 = (AliTPCseed*)array->At(i);
3974     ntracks++;
3975     //
3976     Double_t cradius0 = 40*40;
3977     Double_t cradius1 = 270*270;
3978     Double_t cdist1=8.;
3979     Double_t cdist2=8.;
3980     Double_t cdist3=0.55;
3981     for (Int_t j =i+1;j<nentries;j++){
3982       nall++;
3983       if (sign[j]*sign[i]<1) continue;
3984       if ( (nclusters[i]+nclusters[j])>200) continue;
3985       if ( (nclusters[i]+nclusters[j])<80) continue;
3986       if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
3987       if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
3988       //AliTPCseed * track1 = (AliTPCseed*)array->At(j);  Double_t phase[2][2],radius[2];
3989       Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
3990       if (npoints<1) continue;
3991       // cuts on radius
3992       if (npoints==1){
3993         if (radius[0]<cradius0||radius[0]>cradius1) continue;
3994       }
3995       else{
3996         if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
3997       }
3998       //
3999       Double_t delta1=10000,delta2=10000;
4000       // cuts on the intersection radius
4001       helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4002       if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4003       if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4004       if (npoints==2){
4005         helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4006         if (radius[1]<20&&delta2<1) continue;  //intersection at vertex
4007         if (radius[1]<10&&delta2<3) continue;  //intersection at vertex
4008       }
4009       //
4010       Double_t distance1 = TMath::Min(delta1,delta2);
4011       if (distance1>cdist1) continue;  // cut on DCA linear approximation
4012       //
4013       npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4014       helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4015       if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4016       if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4017       //
4018       if (npoints==2){
4019         helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4020         if (radius[1]<20&&delta2<1) continue;  //intersection at vertex
4021         if (radius[1]<10&&delta2<3) continue;  //intersection at vertex
4022       }
4023       distance1 = TMath::Min(delta1,delta2);
4024       Float_t rkink =0;
4025       if (delta1<delta2){
4026         rkink = TMath::Sqrt(radius[0]);
4027       }
4028       else{
4029         rkink = TMath::Sqrt(radius[1]);
4030       }
4031       if (distance1>cdist2) continue;
4032       //
4033       //
4034       AliTPCseed * track1 = (AliTPCseed*)array->At(j);
4035       //
4036       //
4037       Int_t row0 = GetRowNumber(rkink);
4038       if (row0<10)  continue;
4039       if (row0>150) continue;
4040       //
4041       //
4042       Float_t dens00=-1,dens01=-1;
4043       Float_t dens10=-1,dens11=-1;
4044       //
4045       Int_t found,foundable,shared;
4046       track0->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4047       if (foundable>5) dens00 = Float_t(found)/Float_t(foundable);
4048       track0->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4049       if (foundable>5) dens01 = Float_t(found)/Float_t(foundable);
4050       //
4051       track1->GetClusterStatistic(0,row0-5, found, foundable,shared,kFALSE);
4052       if (foundable>10) dens10 = Float_t(found)/Float_t(foundable);
4053       track1->GetClusterStatistic(row0+5,155, found, foundable,shared,kFALSE);
4054       if (foundable>10) dens11 = Float_t(found)/Float_t(foundable);
4055       //
4056       if (dens00<dens10 && dens01<dens11) continue;
4057       if (dens00>dens10 && dens01>dens11) continue;
4058       if (TMath::Max(dens00,dens10)<0.1)  continue;
4059       if (TMath::Max(dens01,dens11)<0.3)  continue;
4060       //
4061       if (TMath::Min(dens00,dens10)>0.6)  continue;
4062       if (TMath::Min(dens01,dens11)>0.6)  continue;
4063
4064       //
4065       AliTPCseed * ktrack0, *ktrack1;
4066       if (dens00>dens10){
4067         ktrack0 = track0;
4068         ktrack1 = track1;
4069       }
4070       else{
4071         ktrack0 = track1;
4072         ktrack1 = track0;
4073       }
4074       if (TMath::Abs(ktrack0->GetC())>5) continue; // cut on the curvature for mother particle
4075       AliExternalTrackParam paramm(*ktrack0);
4076       AliExternalTrackParam paramd(*ktrack1);
4077       if (row0>60&&ktrack1->GetReference().X()>90.) new (&paramd) AliExternalTrackParam(ktrack1->GetReference());
4078       //
4079       //
4080       kink->SetMother(paramm);
4081       kink->SetDaughter(paramd);
4082       kink->Update();
4083
4084       Float_t x[3] = { kink->fXr[0],kink->fXr[1],kink->fXr[2]};
4085       Int_t index[4];
4086       fParam->Transform0to1(x,index);
4087       fParam->Transform1to2(x,index);
4088       row0 = GetRowNumber(x[0]);
4089
4090       if (kink->fRr<100) continue;
4091       if (kink->fRr>240) continue;
4092       if (kink->fDist2>cdist3) continue;
4093       Float_t dird = kink->fPdr[0]*kink->fXr[0]+kink->fPdr[1]*kink->fXr[1];  // rough direction estimate
4094       if (dird<0) continue;
4095
4096       Float_t dirm = kink->fPm[0]*kink->fXr[0]+kink->fPm[1]*kink->fXr[1];  // rough direction estimate
4097       if (dirm<0) continue;
4098       Float_t mpt = TMath::Sqrt(kink->fPm[0]*kink->fPm[0]+kink->fPm[1]*kink->fPm[1]);
4099       if (mpt<0.2) continue;
4100
4101
4102       Double_t qt   =  TMath::Sin(kink->fAngle[2])*ktrack1->P();
4103       if (qt>0.35) continue;
4104
4105       kink->fLab[0] = CookLabel(ktrack0,0.4,0,row0);
4106       kink->fLab[1] = CookLabel(ktrack1,0.4,row0,160);
4107       if (dens00>dens10){
4108         kink->fTPCdensity[0][0] = dens00;
4109         kink->fTPCdensity[0][1] = dens01;
4110         kink->fTPCdensity[1][0] = dens10;
4111         kink->fTPCdensity[1][1] = dens11;
4112         kink->fIndex[0] = i;
4113         kink->fIndex[1] = j;
4114         kink->fZm[0] = zm[i];
4115         kink->fZm[1] = zm[j];
4116         kink->fFi[0] = fim[i];
4117         kink->fFi[1] = fim[j];
4118       }
4119       else{
4120         kink->fTPCdensity[0][0] = dens10;
4121         kink->fTPCdensity[0][1] = dens11;
4122         kink->fTPCdensity[1][0] = dens00;
4123         kink->fTPCdensity[1][1] = dens01;
4124         kink->fIndex[0] = j;
4125         kink->fIndex[1] = i;
4126         kink->fZm[0] = zm[j];
4127         kink->fZm[1] = zm[i];
4128         kink->fFi[0] = fim[j];
4129         kink->fFi[1] = fim[i];
4130
4131       }
4132       if (kink->GetTPCDensityFactor()<0.8) continue;
4133       if ((2-kink->GetTPCDensityFactor())*kink->fDist2 >0.25) continue;
4134       if (kink->fAngle[2]*ktrack0->P()<0.003) continue; //too small angle
4135       if (kink->fAngle[2]>0.2&&kink->GetTPCDensityFactor()<1.15) continue;
4136       if (kink->fAngle[2]>0.2&&kink->fTPCdensity[0][1]) continue;
4137       if (kink->fAngle[2]<0.02) continue;
4138
4139
4140       //
4141       Int_t drow = Int_t(2./TMath::Tan(0.3+kink->fAngle[2]));  // overlap region defined
4142       kink->fRow0 = row0;
4143       Float_t shapesum =0;
4144       Float_t sum = 0;
4145       for ( Int_t row = row0-drow; row<row0+drow;row++){
4146         if (row<0) continue;
4147         if (row>155) continue;
4148         if (ktrack0->fClusterPointer[row]){
4149           AliTPCTrackerPoint *point =ktrack0->GetTrackPoint(row);
4150           shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4151           sum++;
4152         }
4153         if (ktrack1->fClusterPointer[row]){
4154           AliTPCTrackerPoint *point =ktrack1->GetTrackPoint(row);
4155           shapesum+=point->GetSigmaY()+point->GetSigmaZ();
4156           sum++;
4157         }
4158       }
4159       if (sum<4){
4160         kink->fShapeFactor=-1;
4161       }
4162       else{
4163         kink->fShapeFactor = shapesum/sum;
4164       }
4165
4166       //      esd->AddKink(kink);
4167       kinks->AddLast(kink);
4168       kink = new AliESDkink;
4169       ncandidates++;
4170     }
4171   }
4172   Int_t nkinks = kinks->GetEntriesFast();
4173   Float_t *quality     = new Float_t[nkinks];
4174   Int_t   * indexes = new Int_t[nkinks];
4175   for (Int_t i=0;i<nkinks;i++){
4176     quality[i] =100000;
4177     AliESDkink *kink = (AliESDkink*)kinks->At(i);
4178     if (kink) quality[i] = kink->fDist2*(2.-kink->GetTPCDensityFactor());
4179   }
4180   TMath::Sort(nkinks,quality,indexes,kFALSE);
4181
4182   for (Int_t i=0;i<nkinks;i++){
4183     AliESDkink * kink = (AliESDkink*) kinks->At(indexes[i]);
4184     Int_t index0 = kink->fIndex[0];
4185     Int_t index1 = kink->fIndex[1];
4186     kink->fMultiple[0] = usage[index0];
4187     kink->fMultiple[1] = usage[index1];
4188     if (kink->fMultiple[0]+kink->fMultiple[1]>2) continue;
4189     if (kink->fMultiple[0]+kink->fMultiple[1]>0 && quality[indexes[i]]>0.2) continue;
4190
4191     Int_t index = esd->AddKink(kink);
4192     AliTPCseed * ktrack0 = (AliTPCseed*)array->At(index0);
4193     AliTPCseed * ktrack1 = (AliTPCseed*)array->At(index1);
4194     ktrack0->fKinkIndexes[usage[index0]] = -(index+1);
4195     ktrack1->fKinkIndexes[usage[index1]] =  (index+1);
4196     usage[index0]++;
4197     usage[index1]++;
4198   }
4199   delete [] quality;
4200   delete [] indexes;
4201
4202   delete[] fim;
4203   delete[] zm;
4204   delete [] usage;
4205   delete[] alpha;
4206   delete[] nclusters;
4207   delete[] sign;
4208   delete[] helixes;
4209   kinks->Delete();
4210   delete kinks;
4211
4212   printf("Ncandidates=\t%d\t%d\t%d\n",ncandidates,ntracks,nall);
4213   timer.Print();
4214 }
4215
4216
4217
4218 void AliTPCtrackerMI::UpdateKinkQualityM(AliTPCseed * seed){
4219   //
4220   // update Kink quality information for mother after back propagation
4221   //
4222   if (seed->GetKinkIndex(0)>=0) return;
4223   for (Int_t ikink=0;ikink<3;ikink++){
4224     Int_t index = seed->GetKinkIndex(ikink);
4225     if (index>=0) break;
4226     index = TMath::Abs(index)-1;
4227     AliESDkink * kink = fEvent->GetKink(index);
4228     kink->fTPCdensity2[0][0]=-1;
4229     kink->fTPCdensity2[0][1]=-1;
4230     //
4231     Int_t row0 = kink->fRow0 - Int_t( 2./ (0.05+kink->fAngle[2]));
4232     if (row0<15) row0=15;
4233     //
4234     Int_t row1 = kink->fRow0 + Int_t( 2./ (0.05+kink->fAngle[2]));
4235     if (row1>145) row1=145;
4236     //
4237     Int_t found,foundable,shared;
4238     seed->GetClusterStatistic(0,row0, found, foundable,shared,kFALSE);
4239     if (foundable>5)   kink->fTPCdensity2[0][0] = Float_t(found)/Float_t(foundable);
4240     seed->GetClusterStatistic(row1,155, found, foundable,shared,kFALSE);
4241     if (foundable>5)   kink->fTPCdensity2[0][1] = Float_t(found)/Float_t(foundable);
4242     if (kink->fTPCdensity2[0][1]>0.5) kink->fStatus=-100000;
4243     if (kink->fTPCdensity2[0][0]<0.5) kink->fStatus=-100000;
4244     if (kink->fTPCdensity2[0][0]-kink->fTPCdensity2[0][1]<0.2) kink->fStatus=-100000;
4245
4246     if (kink->fDist2>1) kink->fStatus=-10000;
4247   }
4248
4249 }
4250
4251 Int_t  AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed, Float_t th)
4252 {
4253   //
4254   //
4255   //
4256   for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
4257   //
4258   if (TMath::Abs(seed->GetC())>0.01) return 0;
4259   //
4260
4261   Float_t x[160], y[160], erry[160], z[160], errz[160];
4262   Int_t sec[160];
4263   Float_t xt[160], yt[160], zt[160];
4264   Int_t i1 = 200;
4265   Int_t i2 = 0;
4266   Int_t secm   = -1;
4267   Int_t padm   = -1;
4268   Int_t middle = seed->GetNumberOfClusters()/2;
4269   //
4270   //
4271   // find central sector, get local cooordinates
4272   Int_t count = 0;
4273   for (Int_t i=seed->fFirstPoint;i<=seed->fLastPoint;i++) {
4274     sec[i]= seed->GetClusterSector(i)%18;
4275     x[i]  = GetXrow(i);
4276     if (sec[i]>=0) {
4277       AliTPCclusterMI * cl = seed->fClusterPointer[i];
4278       //      if (cl==0)        cl = GetClusterMI(seed->GetClusterIndex2(i));
4279       if (cl==0) {
4280         sec[i] = -1;
4281         continue;
4282       }
4283       //
4284       //
4285       if (i>i2)  i2 = i;  //last  point with cluster
4286       if (i2<i1) i1 = i;  //first point with cluster
4287       y[i] = cl->GetY();
4288       z[i] = cl->GetZ();
4289       AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
4290       xt[i] = x[i];
4291       yt[i] = point->GetY();
4292       zt[i] = point->GetZ();
4293
4294       if (point->GetX()>0){
4295         erry[i] = point->GetErrY();
4296         errz[i] = point->GetErrZ();
4297       }
4298
4299       count++;
4300       if (count<middle) {
4301         secm = sec[i];  //central sector
4302         padm = i;       //middle point with cluster
4303       }
4304     }
4305   }
4306   //
4307   // rotate position to global coordinate system connected to  sector at last the point
4308   //
4309   for (Int_t i=i1;i<=i2;i++){
4310     //
4311     if (sec[i]<0) continue;
4312     Double_t alpha = (sec[i2]-sec[i])*fSectors->GetAlpha();
4313     Double_t cs = TMath::Cos(alpha);
4314     Double_t sn = TMath::Sin(alpha);
4315     Float_t xx2= x[i]*cs+y[i]*sn;
4316     Float_t yy2= -x[i]*sn+y[i]*cs;
4317     x[i] = xx2;
4318     y[i] = yy2;
4319     //
4320     xx2= xt[i]*cs+yt[i]*sn;
4321     yy2= -xt[i]*sn+yt[i]*cs;
4322     xt[i] = xx2;
4323     yt[i] = yy2;
4324
4325   }
4326   //get "state" vector
4327   Double_t xh[5],xm = x[padm];
4328   xh[0]=yt[i2];
4329   xh[1]=zt[i2];
4330   xh[4]=F1(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
4331   xh[2]=F2(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
4332   xh[3]=F3n(xt[i2],yt[i2],xt[i1],yt[i1],zt[i2],zt[i1],xh[4]);
4333   //
4334   //
4335   for (Int_t i=i1;i<=i2;i++){
4336     Double_t yy,zz;
4337     if (sec[i]<0) continue;
4338     GetProlongation(x[i2], x[i],xh,yy,zz);
4339     if (TMath::Abs(y[i]-yy)>4||TMath::Abs(z[i]-zz)>4){
4340       //Double_t xxh[5];
4341       //xxh[4]=F1old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
4342       //xxh[2]=F2old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
4343       Error("AliTPCtrackerMI::CheckKinkPoint","problem\n");
4344     }
4345     y[i] = y[i] - yy;
4346     z[i] = z[i] - zz;
4347   }
4348   Float_t dyup[160],dydown[160], dzup[160], dzdown[160];
4349   Float_t yup[160], ydown[160],  zup[160],  zdown[160];
4350
4351   AliTPCpolyTrack ptrack1,ptrack2;
4352   //
4353   // derivation up
4354   for (Int_t i=i1;i<=i2;i++){
4355     AliTPCclusterMI * cl = seed->fClusterPointer[i];
4356     if (!cl) continue;
4357     if (cl->GetType()<0) continue;
4358     if (cl->GetType()>10) continue;
4359
4360     if (sec[i]>=0){
4361       ptrack1.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
4362     }
4363     if (ptrack1.GetN()>4.){
4364       ptrack1.UpdateParameters();
4365       Double_t ddy,ddz;
4366       ptrack1.GetFitDerivation(x[i]-xm,ddy,ddz);
4367       Double_t yy,zz;
4368       ptrack1.GetFitPoint(x[i]-xm,yy,zz);
4369
4370       dyup[i] = ddy;
4371       dzup[i] = ddz;
4372       yup[i]  = yy;
4373       zup[i]  = zz;
4374
4375     }
4376     else{
4377       dyup[i]=0.;  //not enough points
4378     }
4379   }
4380   //
4381   // derivation down
4382   for (Int_t i=i2;i>=i1;i--){
4383     AliTPCclusterMI * cl = seed->fClusterPointer[i];
4384     if (!cl) continue;
4385     if (cl->GetType()<0) continue;
4386     if (cl->GetType()>10) continue;
4387     if (sec[i]>=0){
4388       ptrack2.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
4389     }
4390     if (ptrack2.GetN()>4){
4391       ptrack2.UpdateParameters();
4392       Double_t ddy,ddz;
4393       ptrack2.GetFitDerivation(x[i]-xm,ddy,ddz);
4394       Double_t yy,zz;
4395       ptrack2.GetFitPoint(x[i]-xm,yy,zz);
4396
4397       dydown[i] = ddy;
4398       dzdown[i] = ddz;
4399       ydown[i]  = yy;
4400       zdown[i]  = zz;
4401     }
4402     else{
4403       dydown[i]=0.;  //not enough points
4404     }
4405   }
4406   //
4407   //
4408   // find maximal difference of the derivation
4409   for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
4410
4411
4412   for (Int_t i=i1+10;i<i2-10;i++){
4413     if ( (TMath::Abs(dydown[i])<0.00000001)  ||  (TMath::Abs(dyup[i])<0.00000001) ||i<30)continue;
4414     //    printf("%f\t%f\t%f\t%f\t%f\n",x[i],dydown[i],dyup[i],dzdown[i],dzup[i]);
4415     //
4416     Float_t ddy = TMath::Abs(dydown[i]-dyup[i]);
4417     Float_t ddz = TMath::Abs(dzdown[i]-dzup[i]);
4418     if ( (ddy+ddz)> th){
4419       seed->fKinkPoint[0] = i;
4420       seed->fKinkPoint[1] = ddy;
4421       seed->fKinkPoint[2] = ddz;
4422       th = ddy+ddz;
4423     }
4424   }
4425
4426   if (fTreeDebug){
4427     //
4428     //write information to the debug tree
4429     TBranch * br = fTreeDebug->GetBranch("debug");
4430     TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
4431     arr->ExpandCreateFast(i2-i1);
4432     br->SetAddress(&arr);
4433     //
4434     AliTPCclusterMI cldummy;
4435     cldummy.SetQ(0);
4436     AliTPCTrackPoint2 pdummy;
4437     pdummy.GetTPoint().fIsShared = 10;
4438     //
4439     Double_t alpha = sec[i2]*fSectors->GetAlpha();
4440     Double_t cs    = TMath::Cos(alpha);
4441     Double_t sn    = TMath::Sin(alpha);
4442
4443     for (Int_t i=i1;i<i2;i++){
4444       AliTPCTrackPoint2 *trpoint = (AliTPCTrackPoint2*)arr->UncheckedAt(i-i1);
4445       //cluster info
4446       AliTPCclusterMI * cl0 = seed->fClusterPointer[i];
4447       //
4448       AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
4449
4450       if (cl0){
4451         Double_t x = GetXrow(i);
4452         trpoint->GetTPoint() = *point;
4453         trpoint->GetCPoint() = *cl0;
4454         trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
4455         trpoint->fID    = seed->GetUniqueID();
4456         trpoint->fLab   = seed->GetLabel();
4457         //
4458         trpoint->fGX =  cs *x + sn*point->GetY();
4459         trpoint->fGY = -sn *x + cs*point->GetY() ;
4460         trpoint->fGZ = point->GetZ();
4461         //
4462         trpoint->fDY = y[i];
4463         trpoint->fDZ = z[i];
4464         //
4465         trpoint->fDYU = dyup[i];
4466         trpoint->fDZU = dzup[i];
4467         //
4468         trpoint->fDYD = dydown[i];
4469         trpoint->fDZD = dzdown[i];
4470         //
4471         if (TMath::Abs(dyup[i])>0.00000000001 &&TMath::Abs(dydown[i])>0.00000000001){
4472           trpoint->fDDY = dydown[i]-dyup[i];
4473           trpoint->fDDZ = dzdown[i]-dzup[i];
4474         }else{
4475           trpoint->fDDY = 0.;
4476           trpoint->fDDZ = 0.;
4477         }
4478       }
4479       else{
4480         *trpoint = pdummy;
4481         trpoint->GetCPoint()= cldummy;
4482         trpoint->fID = -1;
4483       }
4484       //
4485     }
4486     fTreeDebug->Fill();
4487   }
4488
4489
4490   return 0;
4491
4492 }
4493
4494
4495
4496
4497
4498 AliTPCseed*  AliTPCtrackerMI::ReSeed(AliTPCseed *t)
4499 {
4500   //
4501   // reseed - refit -  track
4502   //
4503   Int_t first = 0;
4504   //  Int_t last  = fSectors->GetNRows()-1;
4505   //
4506   if (fSectors == fOuterSec){
4507     first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows());
4508     //last  =
4509   }
4510   else
4511     first = t->fFirstPoint;
4512   //
4513   AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
4514   FollowBackProlongation(*t,fSectors->GetNRows()-1);
4515   t->Reset(kFALSE);
4516   FollowProlongation(*t,first);
4517   return seed;
4518 }
4519
4520
4521
4522
4523
4524
4525
4526 //_____________________________________________________________________________
4527 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
4528   //-----------------------------------------------------------------
4529   // This function reades track seeds.
4530   //-----------------------------------------------------------------
4531   TDirectory *savedir=gDirectory;
4532
4533   TFile *in=(TFile*)inp;
4534   if (!in->IsOpen()) {
4535      cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
4536      return 1;
4537   }
4538
4539   in->cd();
4540   TTree *seedTree=(TTree*)in->Get("Seeds");
4541   if (!seedTree) {
4542      cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
4543      cerr<<"can't get a tree with track seeds !\n";
4544      return 2;
4545   }
4546   AliTPCtrack *seed=new AliTPCtrack;
4547   seedTree->SetBranchAddress("tracks",&seed);
4548
4549   if (fSeeds==0) fSeeds=new TObjArray(15000);
4550
4551   Int_t n=(Int_t)seedTree->GetEntries();
4552   for (Int_t i=0; i<n; i++) {
4553      seedTree->GetEvent(i);
4554      seed->ResetClusters();
4555      fSeeds->AddLast(new AliTPCseed(*seed/*,seed->GetAlpha()*/));
4556   }
4557
4558   delete seed;
4559   delete seedTree;
4560   savedir->cd();
4561   return 0;
4562 }
4563
4564 Int_t AliTPCtrackerMI::Clusters2Tracks (AliESD *esd)
4565 {
4566   //
4567   if (fSeeds) DeleteSeeds();
4568   fEvent = esd;
4569   Clusters2Tracks();
4570   if (!fSeeds) return 1;
4571   FillESD(fSeeds);
4572   return 0;
4573   //
4574 }
4575
4576
4577 //_____________________________________________________________________________
4578 Int_t AliTPCtrackerMI::Clusters2Tracks() {
4579   //-----------------------------------------------------------------
4580   // This is a track finder.
4581   //-----------------------------------------------------------------
4582   TDirectory *savedir=gDirectory;
4583   TStopwatch timer;
4584
4585   fIteration = 0;
4586   fSeeds = Tracking();
4587
4588   if (fDebug>0){
4589     Info("Clusters2Tracks","Time for tracking: \t");timer.Print();timer.Start();
4590   }
4591   //activate again some tracks
4592   for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
4593     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4594     if (!pt) continue;
4595     Int_t nc=t.GetNumberOfClusters();
4596     if (nc<20) {
4597       delete fSeeds->RemoveAt(i);
4598       continue;
4599     }
4600     if (pt->fRemoval==10) {
4601       if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
4602         pt->Desactivate(10);  // make track again active
4603       else{
4604         pt->Desactivate(20);
4605         delete fSeeds->RemoveAt(i);
4606       }
4607     }
4608   }
4609   //
4610   RemoveUsed2(fSeeds,0.85,0.85,0);
4611   FindKinks(fSeeds,fEvent);
4612   RemoveUsed2(fSeeds,0.5,0.4,30);
4613   //RemoveDouble(fSeeds,0.2,0.6,11);
4614   //  RemoveUsed(fSeeds,0.5,0.5,6);
4615
4616   //
4617   Int_t nseed=fSeeds->GetEntriesFast();
4618   Int_t found = 0;
4619   for (Int_t i=0; i<nseed; i++) {
4620     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4621     if (!pt) continue;
4622     Int_t nc=t.GetNumberOfClusters();
4623     if (nc<15) {
4624       delete fSeeds->RemoveAt(i);
4625       continue;
4626     }
4627     CookLabel(pt,0.1); //For comparison only
4628     //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
4629     if ((pt->IsActive() || (pt->fRemoval==10) )){
4630       found++;
4631       if (fDebug>0) cerr<<found<<'\r';
4632       pt->fLab2 = i;
4633     }
4634     else
4635       delete fSeeds->RemoveAt(i);
4636   }
4637
4638
4639   //RemoveOverlap(fSeeds,0.99,7,kTRUE);
4640   SignShared(fSeeds);
4641   //RemoveUsed(fSeeds,0.9,0.9,6);
4642   //
4643   nseed=fSeeds->GetEntriesFast();
4644   found = 0;
4645   for (Int_t i=0; i<nseed; i++) {
4646     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4647     if (!pt) continue;
4648     Int_t nc=t.GetNumberOfClusters();
4649     if (nc<15) {
4650       delete fSeeds->RemoveAt(i);
4651       continue;
4652     }
4653     t.SetUniqueID(i);
4654     t.CookdEdx(0.02,0.6);
4655     //    CheckKinkPoint(&t,0.05);
4656     //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
4657     if ((pt->IsActive() || (pt->fRemoval==10) )){
4658       found++;
4659       if (fDebug>0){
4660         cerr<<found<<'\r';
4661       }
4662       pt->fLab2 = i;
4663     }
4664     else
4665       delete fSeeds->RemoveAt(i);
4666     //AliTPCseed * seed1 = ReSeed(pt,0.05,0.5,1);
4667     //if (seed1){
4668     //  FollowProlongation(*seed1,0);
4669     //  Int_t n = seed1->GetNumberOfClusters();
4670     //  printf("fP4\t%f\t%f\n",seed1->GetC(),pt->GetC());
4671     //  printf("fN\t%d\t%d\n", seed1->GetNumberOfClusters(),pt->GetNumberOfClusters());
4672     //
4673     //}
4674     //AliTPCseed * seed2 = ReSeed(pt,0.95,0.5,0.05);
4675
4676   }
4677
4678   SortTracks(fSeeds, 1);
4679
4680   /*
4681   fIteration = 1;
4682   PrepareForBackProlongation(fSeeds,5.);
4683   PropagateBack(fSeeds);
4684   printf("Time for back propagation: \t");timer.Print();timer.Start();
4685
4686   fIteration = 2;
4687
4688   PrepareForProlongation(fSeeds,5.);
4689   PropagateForward2(fSeeds);
4690
4691   printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
4692   // RemoveUsed(fSeeds,0.7,0.7,6);
4693   //RemoveOverlap(fSeeds,0.9,7,kTRUE);
4694
4695   nseed=fSeeds->GetEntriesFast();
4696   found = 0;
4697   for (Int_t i=0; i<nseed; i++) {
4698     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
4699     if (!pt) continue;
4700     Int_t nc=t.GetNumberOfClusters();
4701     if (nc<15) {
4702       delete fSeeds->RemoveAt(i);
4703       continue;
4704     }
4705     t.CookdEdx(0.02,0.6);
4706     //    CookLabel(pt,0.1); //For comparison only
4707     //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
4708     if ((pt->IsActive() || (pt->fRemoval==10) )){
4709       cerr<<found++<<'\r';
4710     }
4711     else
4712       delete fSeeds->RemoveAt(i);
4713     pt->fLab2 = i;
4714   }
4715   */
4716
4717   //  fNTracks = found;
4718   if (fDebug>0){
4719     Info("Clusters2Tracks","Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
4720   }
4721   //
4722   //  cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
4723   Info("Clusters2Tracks","Number of found tracks %d",found);
4724   savedir->cd();
4725   //  UnloadClusters();
4726   //
4727   return 0;
4728 }
4729
4730 void AliTPCtrackerMI::Tracking(TObjArray * arr)
4731 {
4732   //
4733   // tracking of the seeds
4734   //
4735
4736   fSectors = fOuterSec;
4737   ParallelTracking(arr,150,63);
4738   fSectors = fOuterSec;
4739   ParallelTracking(arr,63,0);
4740 }
4741
4742 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
4743 {
4744   //
4745   //
4746   //tracking routine
4747   TObjArray * arr = new TObjArray;
4748   //
4749   fSectors = fOuterSec;
4750   TStopwatch timer;
4751   timer.Start();
4752   for (Int_t sec=0;sec<fkNOS;sec++){
4753     if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
4754     if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
4755     if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
4756   }
4757   if (fDebug>0){
4758     Info("Tracking","\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
4759     timer.Print();
4760     timer.Start();
4761   }
4762   Tracking(arr);
4763   if (fDebug>0){
4764     timer.Print();
4765   }
4766
4767   return arr;
4768 }
4769
4770 TObjArray * AliTPCtrackerMI::Tracking()
4771 {
4772   //
4773   //
4774   TStopwatch timer;
4775   timer.Start();
4776   Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
4777
4778   TObjArray * seeds = new TObjArray;
4779   TObjArray * arr=0;
4780
4781   Int_t gap =20;
4782   Float_t cuts[4];
4783   cuts[0] = 0.002;
4784   cuts[1] = 1.5;
4785   cuts[2] = 3.;
4786   cuts[3] = 3.;
4787   Float_t fnumber  = 3.0;
4788   Float_t fdensity = 3.0;
4789
4790   //
4791   //find primaries
4792   cuts[0]=0.0066;
4793   for (Int_t delta = 0; delta<18; delta+=6){
4794     //
4795     cuts[0]=0.0070;
4796     cuts[1] = 1.5;
4797     arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
4798     SumTracks(seeds,arr);
4799     SignClusters(seeds,fnumber,fdensity);
4800     //
4801     for (Int_t i=2;i<6;i+=2){
4802       // seed high pt tracks
4803       cuts[0]=0.0022;
4804       cuts[1]=0.3;
4805       arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
4806       SumTracks(seeds,arr);
4807       SignClusters(seeds,fnumber,fdensity);
4808     }
4809   }
4810   fnumber  = 4;
4811   fdensity = 4.;
4812   //  RemoveUsed(seeds,0.9,0.9,1);
4813   //  UnsignClusters();
4814   //  SignClusters(seeds,fnumber,fdensity);
4815
4816   //find primaries
4817   cuts[0]=0.0077;
4818   for (Int_t delta = 20; delta<120; delta+=10){
4819     //
4820     // seed high pt tracks
4821     cuts[0]=0.0060;
4822     cuts[1]=0.3;
4823     cuts[2]=6.;
4824     arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
4825     SumTracks(seeds,arr);
4826     SignClusters(seeds,fnumber,fdensity);
4827
4828     cuts[0]=0.003;
4829     cuts[1]=0.3;
4830     cuts[2]=6.;
4831     arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
4832     SumTracks(seeds,arr);
4833     SignClusters(seeds,fnumber,fdensity);
4834   }
4835
4836   cuts[0] = 0.01;
4837   cuts[1] = 2.0;
4838   cuts[2] = 3.;
4839   cuts[3] = 2.0;
4840   fnumber  = 2.;
4841   fdensity = 2.;
4842
4843   if (fDebug>0){
4844     Info("Tracking()","\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
4845     timer.Print();
4846     timer.Start();
4847   }
4848   //  RemoveUsed(seeds,0.75,0.75,1);
4849   //UnsignClusters();
4850   //SignClusters(seeds,fnumber,fdensity);
4851
4852   // find secondaries
4853
4854   cuts[0] = 0.3;
4855   cuts[1] = 1.5;
4856   cuts[2] = 3.;
4857   cuts[3] = 1.5;
4858
4859   arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
4860   SumTracks(seeds,arr);
4861   SignClusters(seeds,fnumber,fdensity);
4862   //
4863   arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
4864   SumTracks(seeds,arr);
4865   SignClusters(seeds,fnumber,fdensity);
4866   //
4867   arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
4868   SumTracks(seeds,arr);
4869   SignClusters(seeds,fnumber,fdensity);
4870   //
4871
4872
4873   for (Int_t delta = 3; delta<30; delta+=5){
4874     //
4875     cuts[0] = 0.3;
4876     cuts[1] = 1.5;
4877     cuts[2] = 3.;
4878     cuts[3] = 1.5;
4879     arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4880     SumTracks(seeds,arr);
4881     SignClusters(seeds,fnumber,fdensity);
4882     //
4883     arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
4884     SumTracks(seeds,arr);
4885     SignClusters(seeds,fnumber,fdensity);
4886     //
4887   }
4888   fnumber  = 1;
4889   fdensity = 1;
4890   //
4891   // change cuts
4892   fnumber  = 2.;
4893   fdensity = 2.;
4894   cuts[0]=0.0080;
4895
4896   // find secondaries
4897   for (Int_t delta = 30; delta<70; delta+=10){
4898     //
4899     cuts[0] = 0.3;
4900     cuts[1] = 1.5;
4901     cuts[2] = 3.;
4902     cuts[3] = 1.5;
4903     arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4904     SumTracks(seeds,arr);
4905     SignClusters(seeds,fnumber,fdensity);
4906     //
4907     arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
4908     SumTracks(seeds,arr);
4909     SignClusters(seeds,fnumber,fdensity);
4910   }
4911
4912   if (fDebug>0){
4913     Info("Tracking()","\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
4914     timer.Print();
4915     timer.Start();
4916   }
4917
4918   return seeds;
4919   //
4920
4921 }
4922
4923
4924 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2) const
4925 {
4926   //
4927   //sum tracks to common container
4928   //remove suspicious tracks
4929   Int_t nseed = arr2->GetEntriesFast();
4930   for (Int_t i=0;i<nseed;i++){
4931     AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
4932     if (pt){
4933
4934       // NORMAL ACTIVE TRACK
4935       if (pt->IsActive()){
4936         arr1->AddLast(arr2->RemoveAt(i));
4937         continue;
4938       }
4939       //remove not usable tracks
4940       if (pt->fRemoval!=10){
4941         delete arr2->RemoveAt(i);
4942         continue;
4943       }
4944       // REMOVE VERY SHORT  TRACKS
4945       if (pt->GetNumberOfClusters()<20){
4946         delete arr2->RemoveAt(i);
4947         continue;
4948       }
4949       // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
4950       if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
4951         arr1->AddLast(arr2->RemoveAt(i));
4952       else{
4953         delete arr2->RemoveAt(i);
4954       }
4955     }
4956   }
4957   delete arr2;
4958 }
4959
4960
4961
4962 void  AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
4963 {
4964   //
4965   // try to track in parralel
4966
4967   Int_t nseed=arr->GetEntriesFast();
4968   //prepare seeds for tracking
4969   for (Int_t i=0; i<nseed; i++) {
4970     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
4971     if (!pt) continue;
4972     if (!t.IsActive()) continue;
4973     // follow prolongation to the first layer
4974     if ( (fSectors ==fInnerSec) || (t.fFirstPoint-fParam->GetNRowLow()>rfirst+1) )
4975       FollowProlongation(t, rfirst+1);
4976   }
4977
4978
4979   //
4980   for (Int_t nr=rfirst; nr>=rlast; nr--){
4981     if (nr<fInnerSec->GetNRows())
4982       fSectors = fInnerSec;
4983     else
4984       fSectors = fOuterSec;
4985     // make indexes with the cluster tracks for given
4986
4987     // find nearest cluster
4988     for (Int_t i=0; i<nseed; i++) {
4989       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
4990       if (!pt) continue;
4991       if (nr==80) pt->UpdateReference();
4992       if (!pt->IsActive()) continue;
4993       //      if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
4994       if (pt->fRelativeSector>17) {
4995         continue;
4996       }
4997       UpdateClusters(t,nr);
4998     }
4999     // prolonagate to the nearest cluster - if founded
5000     for (Int_t i=0; i<nseed; i++) {
5001       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
5002       if (!pt) continue;
5003       if (!pt->IsActive()) continue;
5004       // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
5005       if (pt->fRelativeSector>17) {
5006         continue;
5007       }
5008       FollowToNextCluster(*pt,nr);
5009     }
5010   }
5011 }
5012
5013 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac) const
5014 {
5015   //
5016   //
5017   // if we use TPC track itself we have to "update" covariance
5018   //
5019   Int_t nseed= arr->GetEntriesFast();
5020   for (Int_t i=0;i<nseed;i++){
5021     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5022     if (pt) {
5023       pt->Modify(fac);
5024       //
5025       //rotate to current local system at first accepted  point
5026       Int_t index  = pt->GetClusterIndex2(pt->fFirstPoint);
5027       Int_t sec    = (index&0xff000000)>>24;
5028       sec = sec%18;
5029       Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
5030       if (angle1>TMath::Pi())
5031         angle1-=2.*TMath::Pi();
5032       Float_t angle2 = pt->GetAlpha();
5033
5034       if (TMath::Abs(angle1-angle2)>0.001){
5035         pt->Rotate(angle1-angle2);
5036         //angle2 = pt->GetAlpha();
5037         //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
5038         //if (pt->GetAlpha()<0)
5039         //  pt->fRelativeSector+=18;
5040         //sec = pt->fRelativeSector;
5041       }
5042
5043     }
5044
5045   }
5046
5047
5048 }
5049 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac) const
5050 {
5051   //
5052   //
5053   // if we use TPC track itself we have to "update" covariance
5054   //
5055   Int_t nseed= arr->GetEntriesFast();
5056   for (Int_t i=0;i<nseed;i++){
5057     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5058     if (pt) {
5059       pt->Modify(fac);
5060       pt->fFirstPoint = pt->fLastPoint;
5061     }
5062
5063   }
5064
5065
5066 }
5067
5068 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
5069 {
5070   //
5071   // make back propagation
5072   //
5073   Int_t nseed= arr->GetEntriesFast();
5074   for (Int_t i=0;i<nseed;i++){
5075     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5076     if (pt&& pt->GetKinkIndex(0)<=0) {
5077       //AliTPCseed *pt2 = new AliTPCseed(*pt);
5078       fSectors = fInnerSec;
5079       //FollowBackProlongation(*pt,fInnerSec->GetNRows()-1);
5080       //fSectors = fOuterSec;
5081       FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
5082       //if (pt->GetNumberOfClusters()<(pt->fEsd->GetTPCclusters(0)) ){
5083       //        Error("PropagateBack","Not prolonged track %d",pt->GetLabel());
5084       //        FollowBackProlongation(*pt2,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
5085       //}
5086     }
5087     if (pt&& pt->GetKinkIndex(0)>0) {
5088       AliESDkink * kink = fEvent->GetKink(pt->GetKinkIndex(0)-1);
5089       pt->fFirstPoint = kink->fRow0;
5090       fSectors = fInnerSec;
5091       FollowBackProlongation(*pt,fInnerSec->GetNRows()+fOuterSec->GetNRows()-1);
5092     }
5093
5094   }
5095   return 0;
5096 }
5097
5098
5099 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
5100 {
5101   //
5102   // make forward propagation
5103   //
5104   Int_t nseed= arr->GetEntriesFast();
5105   //
5106   for (Int_t i=0;i<nseed;i++){
5107     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
5108     if (pt) {
5109       FollowProlongation(*pt,0);
5110     }
5111   }
5112   return 0;
5113 }
5114
5115
5116 Int_t AliTPCtrackerMI::PropagateForward()
5117 {
5118   //
5119   // propagate track forward
5120   //UnsignClusters();
5121   Int_t nseed = fSeeds->GetEntriesFast();
5122   for (Int_t i=0;i<nseed;i++){
5123     AliTPCseed *pt = (AliTPCseed*)fSeeds->UncheckedAt(i);
5124     if (pt){
5125       AliTPCseed &t = *pt;
5126       Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
5127       if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
5128       if (alpha < 0.            ) alpha += 2.*TMath::Pi();
5129       t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
5130     }
5131   }
5132
5133   fSectors = fOuterSec;
5134   ParallelTracking(fSeeds,fOuterSec->GetNRows()+fInnerSec->GetNRows()-1,fInnerSec->GetNRows());
5135   fSectors = fInnerSec;
5136   ParallelTracking(fSeeds,fInnerSec->GetNRows()-1,0);
5137   //WriteTracks();
5138   return 1;
5139 }
5140
5141
5142
5143
5144
5145
5146 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
5147 {
5148   //
5149   // make back propagation, in between row0 and row1
5150   //
5151
5152   if (pt) {
5153     fSectors = fInnerSec;
5154     Int_t  r1;
5155     //
5156     if (row1<fSectors->GetNRows())
5157       r1 = row1;
5158     else
5159       r1 = fSectors->GetNRows()-1;
5160
5161     if (row0<fSectors->GetNRows()&& r1>0 )
5162       FollowBackProlongation(*pt,r1);
5163     if (row1<=fSectors->GetNRows())
5164       return 0;
5165     //
5166     r1 = row1 - fSectors->GetNRows();
5167     if (r1<=0) return 0;
5168     if (r1>=fOuterSec->GetNRows()) return 0;
5169     fSectors = fOuterSec;
5170     return FollowBackProlongation(*pt,r1);
5171   }
5172   return 0;
5173 }
5174
5175
5176
5177
5178 void  AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
5179 {
5180   //
5181   //
5182   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
5183   //  Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
5184   Float_t padlength =  GetPadPitchLength(row);
5185   //
5186   Float_t sresy = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
5187   Float_t angulary  = seed->GetSnp();
5188   angulary = angulary*angulary/(1-angulary*angulary);
5189   seed->fCurrentSigmaY2 = sd2+padlength*padlength*angulary/12.+sresy*sresy;
5190   //
5191   Float_t sresz = fParam->GetZSigma();
5192   Float_t angularz  = seed->GetTgl();
5193   seed->fCurrentSigmaZ2 = sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz;
5194   /*
5195   Float_t wy = GetSigmaY(seed);
5196   Float_t wz = GetSigmaZ(seed);
5197   wy*=wy;
5198   wz*=wz;
5199   if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
5200     printf("problem\n");
5201   }
5202   */
5203 }
5204
5205
5206 Float_t  AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
5207 {
5208   //
5209   //
5210   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
5211   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
5212   Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
5213   Float_t angular  = seed->GetSnp();
5214   angular = angular*angular/(1-angular*angular);
5215   //  angular*=angular;
5216   //angular  = TMath::Sqrt(angular/(1-angular));
5217   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
5218   return res;
5219 }
5220 Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
5221 {
5222   //
5223   //
5224   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
5225   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
5226   Float_t sres = fParam->GetZSigma();
5227   Float_t angular  = seed->GetTgl();
5228   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
5229   return res;
5230 }
5231
5232
5233 //__________________________________________________________________________
5234 void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
5235   //--------------------------------------------------------------------
5236   //This function "cooks" a track label. If label<0, this track is fake.
5237   //--------------------------------------------------------------------
5238   Int_t noc=t->GetNumberOfClusters();
5239   if (noc<10){
5240     //printf("\nnot founded prolongation\n\n\n");
5241     //t->Dump();
5242     return ;
5243   }
5244   Int_t lb[160];
5245   Int_t mx[160];
5246   AliTPCclusterMI *clusters[160];
5247   //
5248   for (Int_t i=0;i<160;i++) {
5249     clusters[i]=0;
5250     lb[i]=mx[i]=0;
5251   }
5252
5253   Int_t i;
5254   Int_t current=0;
5255   for (i=0; i<160 && current<noc; i++) {
5256
5257      Int_t index=t->GetClusterIndex2(i);
5258      if (index<=0) continue;
5259      if (index&0x8000) continue;
5260      //
5261      //clusters[current]=GetClusterMI(index);
5262      if (t->fClusterPointer[i]){
5263        clusters[current]=t->fClusterPointer[i];
5264        current++;
5265      }
5266   }
5267   noc = current;
5268
5269   Int_t lab=123456789;
5270   for (i=0; i<noc; i++) {
5271     AliTPCclusterMI *c=clusters[i];
5272     if (!c) continue;
5273     lab=TMath::Abs(c->GetLabel(0));
5274     Int_t j;
5275     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
5276     lb[j]=lab;
5277     (mx[j])++;
5278   }
5279
5280   Int_t max=0;
5281   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
5282
5283   for (i=0; i<noc; i++) {
5284     AliTPCclusterMI *c=clusters[i];
5285     if (!c) continue;
5286     if (TMath::Abs(c->GetLabel(1)) == lab ||
5287         TMath::Abs(c->GetLabel(2)) == lab ) max++;
5288   }
5289
5290   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
5291
5292   else {
5293      Int_t tail=Int_t(0.10*noc);
5294      max=0;
5295      Int_t ind=0;
5296      for (i=1; i<=160&&ind<tail; i++) {
5297        //       AliTPCclusterMI *c=clusters[noc-i];
5298        AliTPCclusterMI *c=clusters[i];
5299        if (!c) continue;
5300        if (lab == TMath::Abs(c->GetLabel(0)) ||
5301            lab == TMath::Abs(c->GetLabel(1)) ||
5302            lab == TMath::Abs(c->GetLabel(2))) max++;
5303        ind++;
5304      }
5305      if (max < Int_t(0.5*tail)) lab=-lab;
5306   }
5307
5308   t->SetLabel(lab);
5309
5310   //  delete[] lb;
5311   //delete[] mx;
5312   //delete[] clusters;
5313 }
5314
5315
5316 //__________________________________________________________________________
5317 Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t last) const {
5318   //--------------------------------------------------------------------
5319   //This function "cooks" a track label. If label<0, this track is fake.
5320   //--------------------------------------------------------------------
5321   Int_t noc=t->GetNumberOfClusters();
5322   if (noc<10){
5323     //printf("\nnot founded prolongation\n\n\n");
5324     //t->Dump();
5325     return -1;
5326   }
5327   Int_t lb[160];
5328   Int_t mx[160];
5329   AliTPCclusterMI *clusters[160];
5330   //
5331   for (Int_t i=0;i<160;i++) {
5332     clusters[i]=0;
5333     lb[i]=mx[i]=0;
5334   }
5335
5336   Int_t i;
5337   Int_t current=0;
5338   for (i=0; i<160 && current<noc; i++) {
5339     if (i<first) continue;
5340     if (i>last)  continue;
5341      Int_t index=t->GetClusterIndex2(i);
5342      if (index<=0) continue;
5343      if (index&0x8000) continue;
5344      //
5345      //clusters[current]=GetClusterMI(index);
5346      if (t->fClusterPointer[i]){
5347        clusters[current]=t->fClusterPointer[i];
5348        current++;
5349      }
5350   }
5351   noc = current;
5352   if (noc<5) return -1;
5353   Int_t lab=123456789;
5354   for (i=0; i<noc; i++) {
5355     AliTPCclusterMI *c=clusters[i];
5356     if (!c) continue;
5357     lab=TMath::Abs(c->GetLabel(0));
5358     Int_t j;
5359     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
5360     lb[j]=lab;
5361     (mx[j])++;
5362   }
5363
5364   Int_t max=0;
5365   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
5366
5367   for (i=0; i<noc; i++) {
5368     AliTPCclusterMI *c=clusters[i];
5369     if (!c) continue;
5370     if (TMath::Abs(c->GetLabel(1)) == lab ||
5371         TMath::Abs(c->GetLabel(2)) == lab ) max++;
5372   }
5373
5374   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
5375
5376   else {
5377      Int_t tail=Int_t(0.10*noc);
5378      max=0;
5379      Int_t ind=0;
5380      for (i=1; i<=160&&ind<tail; i++) {
5381        //       AliTPCclusterMI *c=clusters[noc-i];
5382        AliTPCclusterMI *c=clusters[i];
5383        if (!c) continue;
5384        if (lab == TMath::Abs(c->GetLabel(0)) ||
5385            lab == TMath::Abs(c->GetLabel(1)) ||
5386            lab == TMath::Abs(c->GetLabel(2))) max++;
5387        ind++;
5388      }
5389      if (max < Int_t(0.5*tail)) lab=-lab;
5390   }
5391
5392   //  t->SetLabel(lab);
5393   return lab;
5394   //  delete[] lb;
5395   //delete[] mx;
5396   //delete[] clusters;
5397 }
5398
5399
5400 Int_t  AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const
5401 {
5402   //return pad row number for this x
5403   Double_t r;
5404   if (fN < 64){
5405     r=fRow[fN-1].GetX();
5406     if (x > r) return fN;
5407     r=fRow[0].GetX();
5408     if (x < r) return -1;
5409     return Int_t((x-r)/fPadPitchLength + 0.5);}
5410   else{
5411     r=fRow[fN-1].GetX();
5412     if (x > r) return fN;
5413     r=fRow[0].GetX();
5414     if (x < r) return -1;
5415     Double_t r1=fRow[64].GetX();
5416     if(x<r1){
5417       return Int_t((x-r)/f1PadPitchLength + 0.5);}
5418     else{
5419       return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);}
5420   }
5421 }
5422
5423 //_________________________________________________________________________
5424 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
5425   //-----------------------------------------------------------------------
5426   // Setup inner sector
5427   //-----------------------------------------------------------------------
5428   if (f==0) {
5429      fAlpha=par->GetInnerAngle();
5430      fAlphaShift=par->GetInnerAngleShift();
5431      fPadPitchWidth=par->GetInnerPadPitchWidth();
5432      fPadPitchLength=par->GetInnerPadPitchLength();
5433      fN=par->GetNRowLow();
5434      fRow=new AliTPCRow[fN];
5435      for (Int_t i=0; i<fN; i++) {
5436        fRow[i].SetX(par->GetPadRowRadiiLow(i));
5437        fRow[i].fDeadZone =1.5;  //1.5 cm of dead zone
5438      }
5439   } else {
5440      fAlpha=par->GetOuterAngle();
5441      fAlphaShift=par->GetOuterAngleShift();
5442      fPadPitchWidth  = par->GetOuterPadPitchWidth();
5443      fPadPitchLength = par->GetOuter1PadPitchLength();
5444      f1PadPitchLength = par->GetOuter1PadPitchLength();
5445      f2PadPitchLength = par->GetOuter2PadPitchLength();
5446
5447      fN=par->GetNRowUp();
5448      fRow=new AliTPCRow[fN];
5449      for (Int_t i=0; i<fN; i++) {
5450        fRow[i].SetX(par->GetPadRowRadiiUp(i));
5451        fRow[i].fDeadZone =1.5;  // 1.5 cm of dead zone
5452      }
5453   }
5454 }
5455
5456 AliTPCtrackerMI::AliTPCRow::AliTPCRow() {
5457   //
5458   // default constructor
5459   fN=0;
5460   fN1=0;
5461   fN2=0;
5462   fClusters1=0;
5463   fClusters2=0;
5464 }
5465
5466 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
5467   //
5468
5469 }
5470
5471
5472
5473 //_________________________________________________________________________
5474 void
5475 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
5476   //-----------------------------------------------------------------------
5477   // Insert a cluster into this pad row in accordence with its y-coordinate
5478   //-----------------------------------------------------------------------
5479   if (fN==kMaxClusterPerRow) {
5480     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
5481   }
5482   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
5483   Int_t i=Find(c->GetZ());
5484   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
5485   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
5486   fIndex[i]=index; fClusters[i]=c; fN++;
5487 }
5488
5489 void AliTPCtrackerMI::AliTPCRow::ResetClusters() {
5490    //
5491    // reset clusters
5492    fN  = 0;
5493    fN1 = 0;
5494    fN2 = 0;
5495    //delete[] fClusterArray;
5496    if (fClusters1) delete []fClusters1;
5497    if (fClusters2) delete []fClusters2;
5498    //fClusterArray=0;
5499    fClusters1 = 0;
5500    fClusters2 = 0;
5501 }
5502
5503
5504 //___________________________________________________________________
5505 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
5506   //-----------------------------------------------------------------------
5507   // Return the index of the nearest cluster
5508   //-----------------------------------------------------------------------
5509   if (fN==0) return 0;
5510   if (z <= fClusters[0]->GetZ()) return 0;
5511   if (z > fClusters[fN-1]->GetZ()) return fN;
5512   Int_t b=0, e=fN-1, m=(b+e)/2;
5513   for (; b<e; m=(b+e)/2) {
5514     if (z > fClusters[m]->GetZ()) b=m+1;
5515     else e=m;
5516   }
5517   return m;
5518 }
5519
5520
5521
5522 //___________________________________________________________________
5523 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
5524   //-----------------------------------------------------------------------
5525   // Return the index of the nearest cluster in z y
5526   //-----------------------------------------------------------------------
5527   Float_t maxdistance = roady*roady + roadz*roadz;
5528
5529   AliTPCclusterMI *cl =0;
5530   for (Int_t i=Find(z-roadz); i<fN; i++) {
5531       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
5532       if (c->GetZ() > z+roadz) break;
5533       if ( (c->GetY()-y) >  roady ) continue;
5534       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
5535       if (maxdistance>distance) {
5536         maxdistance = distance;
5537         cl=c;
5538       }
5539   }
5540   return cl;
5541 }
5542
5543 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
5544 {
5545   //-----------------------------------------------------------------------
5546   // Return the index of the nearest cluster in z y
5547   //-----------------------------------------------------------------------
5548   Float_t maxdistance = roady*roady + roadz*roadz;
5549   Int_t iz1 = TMath::Max(fFastCluster[Int_t(z-roadz+254.5)]-1,0);
5550   Int_t iz2 = TMath::Min(fFastCluster[Int_t(z+roadz+255.5)]+1,fN);
5551
5552   AliTPCclusterMI *cl =0;
5553   //FindNearest3(y,z,roady,roadz,index);
5554   //  for (Int_t i=Find(z-roadz); i<fN; i++) {
5555   for (Int_t i=iz1; i<iz2; i++) {
5556       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
5557       if (c->GetZ() > z+roadz) break;
5558       if ( c->GetY()-y >  roady ) continue;
5559       if ( y-c->GetY() >  roady ) continue;
5560       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
5561       if (maxdistance>distance) {
5562         maxdistance = distance;
5563         cl=c;
5564         index =i;
5565         //roady = TMath::Sqrt(maxdistance);
5566       }
5567   }
5568   return cl;
5569 }
5570
5571
5572
5573 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
5574 {
5575   //-----------------------------------------------------------------------
5576   // Return the index of the nearest cluster in z y
5577   //-----------------------------------------------------------------------
5578   Float_t maxdistance = roady*roady + roadz*roadz;
5579   //  Int_t iz = Int_t(z+255.);
5580   AliTPCclusterMI *cl =0;
5581   for (Int_t i=Find(z-roadz); i<fN; i++) {
5582     //for (Int_t i=fFastCluster[iz-2]; i<fFastCluster[iz+2]; i++) {
5583       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
5584       if (c->GetZ() > z+roadz) break;
5585       if ( c->GetY()-y >  roady ) continue;
5586       if ( y-c->GetY() >  roady ) continue;
5587       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
5588       if (maxdistance>distance) {
5589         maxdistance = distance;
5590         cl=c;
5591         index =i;
5592         //roady = TMath::Sqrt(maxdistance);
5593       }
5594   }
5595   return cl;
5596 }
5597
5598
5599
5600
5601 AliTPCseed::AliTPCseed():AliTPCtrack(){
5602   //
5603   fRow=0;
5604   fRemoval =0;
5605   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
5606   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
5607   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
5608   fPoints = 0;
5609   fEPoints = 0;
5610   fNFoundable =0;
5611   fNShared  =0;
5612   fRemoval = 0;
5613   fSort =0;
5614   fFirstPoint =0;
5615   fNoCluster =0;
5616   fBSigned = kFALSE;
5617   fSeed1 =-1;
5618   fSeed2 =-1;
5619   fCurrentCluster =0;
5620   fCurrentSigmaY2=0;
5621   fCurrentSigmaZ2=0;
5622 }
5623 AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){
5624   //---------------------
5625   // dummy copy constructor
5626   //-------------------------
5627 }
5628 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
5629   //
5630   //copy constructor
5631   fPoints = 0;
5632   fEPoints = 0;
5633   fNShared  =0;
5634   //  fTrackPoints =0;
5635   fRemoval =0;
5636   fSort =0;
5637   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=t.GetKinkIndex(i);
5638
5639   for (Int_t i=0;i<160;i++) {
5640     fClusterPointer[i] = 0;
5641     Int_t index = t.GetClusterIndex(i);
5642     if (index>=-1){
5643       SetClusterIndex2(i,index);
5644     }
5645     else{
5646       SetClusterIndex2(i,-3);
5647     }
5648   }
5649   fFirstPoint =0;
5650   fNoCluster =0;
5651   fBSigned = kFALSE;
5652   fSeed1 =-1;
5653   fSeed2 =-1;
5654   fCurrentCluster =0;
5655   fCurrentSigmaY2=0;
5656   fCurrentSigmaZ2=0;
5657 }
5658 /*
5659 AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
5660   //
5661   //copy constructor
5662   fRow=0;
5663   for (Int_t i=0;i<160;i++) {
5664     fClusterPointer[i] = 0;
5665     Int_t index = t.GetClusterIndex(i);
5666     SetClusterIndex2(i,index);
5667   }
5668   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
5669   fPoints = 0;
5670   fEPoints = 0;
5671   fNFoundable =0;
5672   fNShared  =0;
5673   //  fTrackPoints =0;
5674   fRemoval =0;
5675   fSort = 0;
5676   fFirstPoint =0;
5677   fNoCluster =0;
5678   fBSigned = kFALSE;
5679   fSeed1 =-1;
5680   fSeed2 =-1;
5681   fCurrentCluster =0;
5682   fCurrentSigmaY2=0;
5683   fCurrentSigmaZ2=0;
5684
5685 }
5686 */
5687
5688 AliTPCseed::AliTPCseed(UInt_t index,  const Double_t xx[5], const Double_t cc[15],
5689                                         Double_t xr, Double_t alpha):
5690   AliTPCtrack(index, xx, cc, xr, alpha) {
5691    //
5692   //
5693   //constructor
5694   fRow =0;
5695   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
5696   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
5697   for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
5698   fPoints = 0;
5699   fEPoints = 0;
5700   fNFoundable =0;
5701   fNShared  = 0;
5702   //  fTrackPoints =0;
5703   fRemoval =0;
5704   fSort =0;
5705   fFirstPoint =0;
5706   //  fHelixIn = new TClonesArray("AliHelix",0);
5707   //fHelixOut = new TClonesArray("AliHelix",0);
5708   fNoCluster =0;
5709   fBSigned = kFALSE;
5710   fSeed1 =-1;
5711   fSeed2 =-1;
5712   fCurrentCluster =0;
5713   fCurrentSigmaY2=0;
5714   fCurrentSigmaZ2=0;
5715 }
5716
5717 AliTPCseed::~AliTPCseed(){
5718   //
5719   // destructor
5720   if (fPoints) delete fPoints;
5721   fPoints =0;
5722   if (fEPoints) delete fEPoints;
5723   fEPoints = 0;
5724   fNoCluster =0;
5725 }
5726
5727 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
5728 {
5729   //
5730   //
5731   return &fTrackPoints[i];
5732 }
5733
5734 void AliTPCseed::RebuildSeed()
5735 {
5736   //
5737   // rebuild seed to be ready for storing
5738   AliTPCclusterMI cldummy;
5739   cldummy.SetQ(0);
5740   AliTPCTrackPoint pdummy;
5741   pdummy.GetTPoint().fIsShared = 10;
5742   for (Int_t i=0;i<160;i++){
5743     AliTPCclusterMI * cl0 = fClusterPointer[i];
5744     AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
5745     if (cl0){
5746       trpoint->GetTPoint() = *(GetTrackPoint(i));
5747       trpoint->GetCPoint() = *cl0;
5748       trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
5749     }
5750     else{
5751       *trpoint = pdummy;
5752       trpoint->GetCPoint()= cldummy;
5753     }
5754
5755   }
5756
5757 }
5758
5759
5760 Double_t AliTPCseed::GetDensityFirst(Int_t n)
5761 {
5762   //
5763   //
5764   // return cluster for n rows bellow first point
5765   Int_t nfoundable = 1;
5766   Int_t nfound      = 1;
5767   for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
5768     Int_t index = GetClusterIndex2(i);
5769     if (index!=-1) nfoundable++;
5770     if (index>0) nfound++;
5771   }
5772   if (nfoundable<n) return 0;
5773   return Double_t(nfound)/Double_t(nfoundable);
5774
5775 }
5776
5777
5778 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
5779 {
5780   // get cluster stat.  on given region
5781   //
5782   found       = 0;
5783   foundable   = 0;
5784   shared      =0;
5785   for (Int_t i=first;i<last; i++){
5786     Int_t index = GetClusterIndex2(i);
5787     if (index!=-1) foundable++;
5788     if (fClusterPointer[i]) {
5789       found++;
5790     }
5791     else
5792       continue;
5793
5794     if (fClusterPointer[i]->IsUsed(10)) {
5795       shared++;
5796       continue;
5797     }
5798     if (!plus2) continue; //take also neighborhoud
5799     //
5800     if ( (i>0) && fClusterPointer[i-1]){
5801       if (fClusterPointer[i-1]->IsUsed(10)) {
5802         shared++;
5803         continue;
5804       }
5805     }
5806     if ( fClusterPointer[i+1]){
5807       if (fClusterPointer[i+1]->IsUsed(10)) {
5808         shared++;
5809         continue;
5810       }
5811     }
5812
5813   }
5814   //if (shared>found){
5815     //Error("AliTPCseed::GetClusterStatistic","problem\n");
5816   //}
5817 }
5818
5819 //_____________________________________________________________________________
5820 void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
5821   //-----------------------------------------------------------------
5822   // This funtion calculates dE/dX within the "low" and "up" cuts.
5823   //-----------------------------------------------------------------
5824
5825   Float_t amp[200];
5826   Float_t angular[200];
5827   Float_t weight[200];
5828   Int_t index[200];
5829   //Int_t nc = 0;
5830   //  TClonesArray & arr = *fPoints;
5831   Float_t meanlog = 100.;
5832
5833   Float_t mean[4]  = {0,0,0,0};
5834   Float_t sigma[4] = {1000,1000,1000,1000};
5835   Int_t nc[4]      = {0,0,0,0};
5836   Float_t norm[4]    = {1000,1000,1000,1000};
5837   //
5838   //
5839   fNShared =0;
5840
5841   for (Int_t of =0; of<4; of++){
5842     for (Int_t i=of+i1;i<i2;i+=4)
5843       {
5844         Int_t index = fIndex[i];
5845         if (index<0||index&0x8000) continue;
5846
5847         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
5848         AliTPCTrackerPoint * point = GetTrackPoint(i);
5849         //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
5850         //AliTPCTrackerPoint * pointp = 0;
5851         //if (i<159) pointp = GetTrackPoint(i+1);
5852
5853         if (point==0) continue;
5854         AliTPCclusterMI * cl = fClusterPointer[i];
5855         if (cl==0) continue;
5856         if (onlyused && (!cl->IsUsed(10))) continue;
5857         if (cl->IsUsed(11)) {
5858           fNShared++;
5859           continue;
5860         }
5861         Int_t   type   = cl->GetType();
5862         //if (point->fIsShared){
5863         //  fNShared++;
5864         //  continue;
5865         //}
5866         //if (pointm)
5867         //  if (pointm->fIsShared) continue;
5868         //if (pointp)
5869         //  if (pointp->fIsShared) continue;
5870
5871         if (type<0) continue;
5872         //if (type>10) continue;
5873         //if (point->GetErrY()==0) continue;
5874         //if (point->GetErrZ()==0) continue;
5875
5876         //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
5877         //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
5878         //if ((ddy*ddy+ddz*ddz)>10) continue;
5879
5880
5881         //      if (point->GetCPoint().GetMax()<5) continue;
5882         if (cl->GetMax()<5) continue;
5883         Float_t angley = point->GetAngleY();
5884         Float_t anglez = point->GetAngleZ();
5885
5886         Float_t rsigmay2 =  point->GetSigmaY();
5887         Float_t rsigmaz2 =  point->GetSigmaZ();
5888         /*
5889         Float_t ns = 1.;
5890         if (pointm){
5891           rsigmay +=  pointm->GetTPoint().GetSigmaY();
5892           rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
5893           ns+=1.;
5894         }
5895         if (pointp){
5896           rsigmay +=  pointp->GetTPoint().GetSigmaY();
5897           rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
5898           ns+=1.;
5899         }
5900         rsigmay/=ns;
5901         rsigmaz/=ns;
5902         */
5903
5904         Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
5905
5906         Float_t ampc   = 0;     // normalization to the number of electrons
5907         if (i>64){
5908           //      ampc = 1.*point->GetCPoint().GetMax();
5909           ampc = 1.*cl->GetMax();
5910           //ampc = 1.*point->GetCPoint().GetQ();
5911           //      AliTPCClusterPoint & p = point->GetCPoint();
5912           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
5913           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
5914           //Float_t dz =
5915           //  TMath::Abs( Int_t(iz) - iz + 0.5);
5916           //ampc *= 1.15*(1-0.3*dy);
5917           //ampc *= 1.15*(1-0.3*dz);
5918           //      Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
5919           //ampc               *=zfactor;
5920         }
5921         else{
5922           //ampc = 1.0*point->GetCPoint().GetMax();
5923           ampc = 1.0*cl->GetMax();
5924           //ampc = 1.0*point->GetCPoint().GetQ();
5925           //AliTPCClusterPoint & p = point->GetCPoint();
5926           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
5927           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
5928           //Float_t dz =
5929           //  TMath::Abs( Int_t(iz) - iz + 0.5);
5930
5931           //ampc *= 1.15*(1-0.3*dy);
5932           //ampc *= 1.15*(1-0.3*dz);
5933           //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
5934           //ampc               *=zfactor;
5935
5936         }
5937         ampc *= 2.0;     // put mean value to channel 50
5938         //ampc *= 0.58;     // put mean value to channel 50
5939         Float_t w      =  1.;
5940         //      if (type>0)  w =  1./(type/2.-0.5);
5941         //      Float_t z = TMath::Abs(cl->GetZ());
5942         if (i<64) {
5943           ampc /= 0.6;
5944           //ampc /= (1+0.0008*z);
5945         } else
5946           if (i>128){
5947             ampc /=1.5;
5948             //ampc /= (1+0.0008*z);
5949           }else{
5950             //ampc /= (1+0.0008*z);
5951           }
5952
5953         if (type<0) {  //amp at the border - lower weight
5954           // w*= 2.;
5955
5956           continue;
5957         }
5958         if (rsigma>1.5) ampc/=1.3;  // if big backround
5959         amp[nc[of]]        = ampc;
5960         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
5961         weight[nc[of]]     = w;
5962         nc[of]++;
5963       }
5964
5965     TMath::Sort(nc[of],amp,index,kFALSE);
5966     Float_t sumamp=0;
5967     Float_t sumamp2=0;
5968     Float_t sumw=0;
5969     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
5970     meanlog = 50;
5971     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
5972       Float_t ampl      = amp[index[i]]/angular[index[i]];
5973       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
5974       //
5975       sumw    += weight[index[i]];
5976       sumamp  += weight[index[i]]*ampl;
5977       sumamp2 += weight[index[i]]*ampl*ampl;
5978       norm[of]    += angular[index[i]]*weight[index[i]];
5979     }
5980     if (sumw<1){
5981       SetdEdx(0);
5982     }
5983     else {
5984       norm[of] /= sumw;
5985       mean[of]  = sumamp/sumw;
5986       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
5987       if (sigma[of]>0.1)
5988         sigma[of] = TMath::Sqrt(sigma[of]);
5989       else
5990         sigma[of] = 1000;
5991
5992     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
5993     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
5994     //mean *=(1-0.1*(norm-1.));
5995     }
5996   }
5997
5998   Float_t dedx =0;
5999   fSdEdx =0;
6000   fMAngular =0;
6001   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
6002   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
6003
6004
6005   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
6006   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
6007
6008   Int_t norm2 = 0;
6009   Int_t norm3 = 0;
6010   for (Int_t i =0;i<4;i++){
6011     if (nc[i]>2&&nc[i]<1000){
6012       dedx      += mean[i] *nc[i];
6013       fSdEdx    += sigma[i]*(nc[i]-2);
6014       fMAngular += norm[i] *nc[i];
6015       norm2     += nc[i];
6016       norm3     += nc[i]-2;
6017     }
6018     fDEDX[i]  = mean[i];
6019     fSDEDX[i] = sigma[i];
6020     fNCDEDX[i]= nc[i];
6021   }
6022
6023   if (norm3>0){
6024     dedx   /=norm2;
6025     fSdEdx /=norm3;
6026     fMAngular/=norm2;
6027   }
6028   else{
6029     SetdEdx(0);
6030     return;
6031   }
6032   //  Float_t dedx1 =dedx;
6033   /*
6034   dedx =0;
6035   for (Int_t i =0;i<4;i++){
6036     if (nc[i]>2&&nc[i]<1000){
6037       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
6038       dedx      += mean[i] *nc[i];
6039     }
6040     fDEDX[i]  = mean[i];
6041   }
6042   dedx /= norm2;
6043   */
6044
6045
6046   SetdEdx(dedx);
6047
6048   //mi deDX
6049
6050
6051
6052   //Very rough PID
6053   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
6054
6055   if (p<0.6) {
6056     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(AliPID::ParticleMass(AliPID::kPion)); return;}
6057     if (dedx < 39.+ 12./p/p) { SetMass(AliPID::ParticleMass(AliPID::kKaon)); return;}
6058     SetMass(AliPID::ParticleMass(AliPID::kProton)); return;
6059   }
6060
6061   if (p<1.2) {
6062     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(AliPID::ParticleMass(AliPID::kPion)); return;}
6063     SetMass(AliPID::ParticleMass(AliPID::kProton)); return;
6064   }
6065
6066   SetMass(AliPID::ParticleMass(AliPID::kPion)); return;
6067
6068 }
6069
6070
6071
6072 /*
6073
6074
6075
6076 void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
6077   //-----------------------------------------------------------------
6078   // This funtion calculates dE/dX within the "low" and "up" cuts.
6079   //-----------------------------------------------------------------
6080
6081   Float_t amp[200];
6082   Float_t angular[200];
6083   Float_t weight[200];
6084   Int_t index[200];
6085   Bool_t inlimit[200];
6086   for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
6087   for (Int_t i=0;i<200;i++) amp[i]=10000;
6088   for (Int_t i=0;i<200;i++) angular[i]= 1;;
6089
6090
6091   //
6092   Float_t meanlog = 100.;
6093   Int_t indexde[4]={0,64,128,160};
6094
6095   Float_t amean     =0;
6096   Float_t asigma    =0;
6097   Float_t anc       =0;
6098   Float_t anorm     =0;
6099
6100   Float_t mean[4]  = {0,0,0,0};
6101   Float_t sigma[4] = {1000,1000,1000,1000};
6102   Int_t nc[4]      = {0,0,0,0};
6103   Float_t norm[4]    = {1000,1000,1000,1000};
6104   //
6105   //
6106   fNShared =0;
6107
6108   //  for (Int_t of =0; of<3; of++){
6109   //  for (Int_t i=indexde[of];i<indexde[of+1];i++)
6110   for (Int_t i =0; i<160;i++)
6111     {
6112         AliTPCTrackPoint * point = GetTrackPoint(i);
6113         if (point==0) continue;
6114         if (point->fIsShared){
6115           fNShared++;
6116           continue;
6117         }
6118         Int_t   type   = point->GetCPoint().GetType();
6119         if (type<0) continue;
6120         if (point->GetCPoint().GetMax()<5) continue;
6121         Float_t angley = point->GetTPoint().GetAngleY();
6122         Float_t anglez = point->GetTPoint().GetAngleZ();
6123         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
6124         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
6125         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
6126
6127         Float_t ampc   = 0;     // normalization to the number of electrons
6128         if (i>64){
6129           ampc =  point->GetCPoint().GetMax();
6130         }
6131         else{
6132           ampc = point->GetCPoint().GetMax();
6133         }
6134         ampc *= 2.0;     // put mean value to channel 50
6135         //      ampc *= 0.565;     // put mean value to channel 50
6136
6137         Float_t w      =  1.;
6138         Float_t z = TMath::Abs(point->GetCPoint().GetZ());
6139         if (i<64) {
6140           ampc /= 0.63;
6141         } else
6142           if (i>128){
6143             ampc /=1.51;
6144           }
6145         if (type<0) {  //amp at the border - lower weight
6146           continue;
6147         }
6148         if (rsigma>1.5) ampc/=1.3;  // if big backround
6149         angular[i]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
6150         amp[i]        = ampc/angular[i];
6151         weight[i]     = w;
6152         anc++;
6153     }
6154
6155   TMath::Sort(159,amp,index,kFALSE);
6156   for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
6157     inlimit[index[i]] = kTRUE;  // take all clusters
6158   }
6159
6160   //  meanlog = amp[index[Int_t(anc*0.3)]];
6161   meanlog =10000.;
6162   for (Int_t of =0; of<3; of++){
6163     Float_t sumamp=0;
6164     Float_t sumamp2=0;
6165     Float_t sumw=0;
6166    for (Int_t i=indexde[of];i<indexde[of+1];i++)
6167       {
6168         if (inlimit[i]==kFALSE) continue;
6169         Float_t ampl      = amp[i];
6170         ///angular[i];
6171         ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
6172         //
6173         sumw    += weight[i];
6174         sumamp  += weight[i]*ampl;
6175         sumamp2 += weight[i]*ampl*ampl;
6176         norm[of]    += angular[i]*weight[i];
6177         nc[of]++;
6178       }
6179    if (sumw<1){
6180      SetdEdx(0);
6181    }
6182    else {
6183      norm[of] /= sumw;
6184      mean[of]  = sumamp/sumw;
6185      sigma[of] = sumamp2/sumw-mean[of]*mean[of];
6186      if (sigma[of]>0.1)
6187        sigma[of] = TMath::Sqrt(sigma[of]);
6188      else
6189        sigma[of] = 1000;
6190      mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
6191    }
6192   }
6193
6194   Float_t dedx =0;
6195   fSdEdx =0;
6196   fMAngular =0;
6197   //
6198   Int_t norm2 = 0;
6199   Int_t norm3 = 0;
6200   Float_t www[3] = {12.,14.,17.};
6201   //Float_t www[3] = {1.,1.,1.};
6202
6203   for (Int_t i =0;i<3;i++){
6204     if (nc[i]>2&&nc[i]<1000){
6205       dedx      += mean[i] *nc[i]*www[i]/sigma[i];
6206       fSdEdx    += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
6207       fMAngular += norm[i] *nc[i];
6208       norm2     += nc[i]*www[i]/sigma[i];
6209       norm3     += (nc[i]-2)*www[i]/sigma[i];
6210     }
6211     fDEDX[i]  = mean[i];
6212     fSDEDX[i] = sigma[i];
6213     fNCDEDX[i]= nc[i];
6214   }
6215
6216   if (norm3>0){
6217     dedx   /=norm2;
6218     fSdEdx /=norm3;
6219     fMAngular/=norm2;
6220   }
6221   else{
6222     SetdEdx(0);
6223     return;
6224   }
6225   //  Float_t dedx1 =dedx;
6226
6227   dedx =0;
6228   Float_t norm4 = 0;
6229   for (Int_t i =0;i<3;i++){
6230     if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
6231       //mean[i]   = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
6232       dedx      += mean[i] *(nc[i])/(sigma[i]);
6233       norm4     += (nc[i])/(sigma[i]);
6234     }
6235     fDEDX[i]  = mean[i];
6236   }
6237   if (norm4>0) dedx /= norm4;
6238
6239
6240
6241   SetdEdx(dedx);
6242
6243   //mi deDX
6244
6245 }
6246
6247 */