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 #include "Riostream.h"
  34 #include <TClonesArray.h>
  35 #include <TFile.h>
  36 #include <TObjArray.h>
  37 #include <TTree.h>
  38
  39 #include "AliComplexCluster.h"
  40 #include "AliESD.h"
  41 #include "AliESDkink.h"
  42 #include "AliHelix.h"
  43 #include "AliRunLoader.h"
  44 #include "AliTPCClustersRow.h"
  45 #include "AliTPCParam.h"
  46 #include "AliTPCReconstructor.h"
  47 #include "AliTPCclusterMI.h"
  48 #include "AliTPCpolyTrack.h"
  49 #include "AliTPCreco.h"
  50 #include "AliTPCtrackerMI.h"
  51 #include "TStopwatch.h"
  52 #include "TTreeStream.h"
  53 //
  54
  55 ClassImp(AliTPCseed)
  56 ClassImp(AliTPCtrackerMI)
  57
  58
  59 class AliTPCFastMath {
  60 public:
  61   AliTPCFastMath();
  62   static Double_t FastAsin(Double_t x);
  63  private:
  64   static Double_t fgFastAsin[20000];  //lookup table for fast asin computation
  65 };
  66
  67 Double_t AliTPCFastMath::fgFastAsin[20000];
  68 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
  69
  70 AliTPCFastMath::AliTPCFastMath(){
  71   //
  72   // initialized lookup table;
  73   for (Int_t i=0;i<10000;i++){
  74     fgFastAsin[2*i] = TMath::ASin(i/10000.);
  75     fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
  76   }
  77 }
  78
  79 Double_t AliTPCFastMath::FastAsin(Double_t x){
  80   //
  81   // return asin using lookup table
  82   if (x>0){
  83     Int_t index = int(x*10000);
  84     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
  85   }
  86   x*=-1;
  87   Int_t index = int(x*10000);
  88   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
  89 }
  90
  91
  92
  93
  94 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
  95   //
  96   //update track information using current cluster - track->fCurrentCluster
  97
  98
  99   AliTPCclusterMI* c =track->fCurrentCluster;
 100   if (accept>0) track->fCurrentClusterIndex1 |=0x8000;  //sign not accepted clusters
 101
 102   UInt_t i = track->fCurrentClusterIndex1;
 103
 104   Int_t sec=(i&0xff000000)>>24;
 105   //Int_t row = (i&0x00ff0000)>>16;
 106   track->fRow=(i&0x00ff0000)>>16;
 107   track->fSector = sec;
 108   //  Int_t index = i&0xFFFF;
 109   if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
 110   track->SetClusterIndex2(track->fRow, i);
 111   //track->fFirstPoint = row;
 112   //if ( track->fLastPoint<row) track->fLastPoint =row;
 113   //  if (track->fRow<0 || track->fRow>160) {
 114   //  printf("problem\n");
 115   //}
 116   if (track->fFirstPoint>track->fRow)
 117     track->fFirstPoint = track->fRow;
 118   if (track->fLastPoint<track->fRow)
 119     track->fLastPoint  = track->fRow;
 120
 121
 122   track->fClusterPointer[track->fRow] = c;
 123   //
 124
 125   Float_t angle2 = track->GetSnp()*track->GetSnp();
 126   angle2 = TMath::Sqrt(angle2/(1-angle2));
 127   //
 128   //SET NEW Track Point
 129   //
 130   //  if (debug)
 131   {
 132     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->fRow));
 133     //
 134     point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
 135     point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
 136     point.SetErrY(sqrt(track->fErrorY2));
 137     point.SetErrZ(sqrt(track->fErrorZ2));
 138     //
 139     point.SetX(track->GetX());
 140     point.SetY(track->GetY());
 141     point.SetZ(track->GetZ());
 142     point.SetAngleY(angle2);
 143     point.SetAngleZ(track->GetTgl());
 144     if (point.fIsShared){
 145       track->fErrorY2 *= 4;
 146       track->fErrorZ2 *= 4;
 147     }
 148   }
 149
 150   Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
 151   //
 152   track->fErrorY2 *= 1.3;
 153   track->fErrorY2 += 0.01;
 154   track->fErrorZ2 *= 1.3;
 155   track->fErrorZ2 += 0.005;
 156     //}
 157   if (accept>0) return 0;
 158   if (track->GetNumberOfClusters()%20==0){
 159     //    if (track->fHelixIn){
 160     //  TClonesArray & larr = *(track->fHelixIn);
 161     //  Int_t ihelix = larr.GetEntriesFast();
 162     //  new(larr[ihelix]) AliHelix(*track) ;
 163     //}
 164   }
 165   track->fNoCluster =0;
 166   return track->Update(c,chi2,i);
 167 }
 168
 169
 170
 171 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
 172                                       Float_t cory, Float_t corz)
 173 {
 174   //
 175   // decide according desired precision to accept given
 176   // cluster for tracking
 177   Double_t sy2=ErrY2(seed,cluster)*cory;
 178   Double_t sz2=ErrZ2(seed,cluster)*corz;
 179   //sy2=ErrY2(seed,cluster)*cory;
 180   //sz2=ErrZ2(seed,cluster)*cory;
 181
 182   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
 183   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
 184
 185   Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
 186     (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
 187   Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
 188     (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
 189
 190   Double_t rdistance2  = rdistancey2+rdistancez2;
 191   //Int_t  accept =0;
 192
 193   if (rdistance2>16) return 3;
 194
 195
 196   if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
 197     return 2;  //suspisiouce - will be changed
 198
 199   if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
 200     // strict cut on overlaped cluster
 201     return  2;  //suspisiouce - will be changed
 202
 203   if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
 204        && cluster->GetType()<0){
 205     seed->fNFoundable--;
 206     return 2;
 207   }
 208   return 0;
 209 }
 210
 211
 212
 213
 214 //_____________________________________________________________________________
 215 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
 216 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 217 {
 218   //---------------------------------------------------------------------
 219   // The main TPC tracker constructor
 220   //---------------------------------------------------------------------
 221   fInnerSec=new AliTPCSector[fkNIS];
 222   fOuterSec=new AliTPCSector[fkNOS];
 223
 224   Int_t i;
 225   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
 226   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
 227
 228   fN=0;  fSectors=0;
 229
 230   fSeeds=0;
 231   fNtracks = 0;
 232   fParam = par;
 233   Int_t nrowlow = par->GetNRowLow();
 234   Int_t nrowup = par->GetNRowUp();
 235
 236
 237   for (Int_t i=0;i<nrowlow;i++){
 238     fXRow[i]     = par->GetPadRowRadiiLow(i);
 239     fPadLength[i]= par->GetPadPitchLength(0,i);
 240     fYMax[i]     = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
 241   }
 242
 243
 244   for (Int_t i=0;i<nrowup;i++){
 245     fXRow[i+nrowlow]      = par->GetPadRowRadiiUp(i);
 246     fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
 247     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
 248   }
 249   fSeeds=0;
 250   //
 251   fInput    = 0;
 252   fOutput   = 0;
 253   fSeedTree = 0;
 254   fTreeDebug =0;
 255   fNewIO     =0;
 256   fDebug     =0;
 257   fEvent     =0;
 258 }
 259 //________________________________________________________________________
 260 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
 261   AliTracker(t),
 262   fkNIS(t.fkNIS),
 263   fkNOS(t.fkNOS)
 264 {
 265   //------------------------------------
 266   // dummy copy constructor
 267   //------------------------------------------------------------------
 268 }
 269 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
 270   //------------------------------
 271   // dummy
 272   //--------------------------------------------------------------
 273   return *this;
 274 }
 275 //_____________________________________________________________________________
 276 AliTPCtrackerMI::~AliTPCtrackerMI() {
 277   //------------------------------------------------------------------
 278   // TPC tracker destructor
 279   //------------------------------------------------------------------
 280   delete[] fInnerSec;
 281   delete[] fOuterSec;
 282   if (fSeeds) {
 283     fSeeds->Delete();
 284     delete fSeeds;
 285   }
 286 }
 287
 288 void AliTPCtrackerMI::SetIO()
 289 {
 290   //
 291   fNewIO   =  kTRUE;
 292   fInput   =  AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::GetDefaultEventFolderName());
 293
 294   fOutput  =  AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::GetDefaultEventFolderName());
 295   if (fOutput){
 296     AliTPCtrack *iotrack= new AliTPCtrack;
 297     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 298     delete iotrack;
 299   }
 300 }
 301
 302
 303 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
 304 {
 305
 306   // set input
 307   fNewIO = kFALSE;
 308   fInput    = 0;
 309   fOutput   = 0;
 310   fSeedTree = 0;
 311   fTreeDebug =0;
 312   fInput = input;
 313   if (input==0){
 314     return;
 315   }
 316   //set output
 317   fOutput = output;
 318   if (output){
 319     AliTPCtrack *iotrack= new AliTPCtrack;
 320     //    iotrack->fHelixIn   = new TClonesArray("AliHelix");
 321     //iotrack->fHelixOut  = new TClonesArray("AliHelix");
 322     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 323     delete iotrack;
 324   }
 325   if (output && (fDebug&2)){
 326     //write the full seed information if specified in debug mode
 327     //
 328     fSeedTree =  new TTree("Seeds","Seeds");
 329     AliTPCseed * vseed = new AliTPCseed;
 330     //
 331     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
 332     arrtr->ExpandCreateFast(160);
 333     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
 334     //
 335     vseed->fPoints = arrtr;
 336     vseed->fEPoints = arre;
 337     //    vseed->fClusterPoints = arrcl;
 338     fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
 339     delete arrtr;
 340     delete arre;
 341     fTreeDebug = new TTree("trackDebug","trackDebug");
 342     TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
 343     fTreeDebug->Branch("debug",&arrd,32000,99);
 344   }
 345
 346
 347   //set ESD event
 348   fEvent  = event;
 349 }
 350
 351 void AliTPCtrackerMI::FillESD(TObjArray* arr)
 352 {
 353   //
 354   //
 355   //fill esds using updated tracks
 356   if (fEvent){
 357     // write tracks to the event
 358     // store index of the track
 359     Int_t nseed=arr->GetEntriesFast();
 360     //FindKinks(arr,fEvent);
 361     for (Int_t i=0; i<nseed; i++) {
 362       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
 363       if (!pt) continue;
 364       pt->UpdatePoints();
 365       //      pt->PropagateTo(fParam->GetInnerRadiusLow());
 366       if (pt->GetKinkIndex(0)<=0){  //don't propagate daughter tracks
 367         pt->PropagateTo(fParam->GetInnerRadiusLow());
 368       }
 369
 370       if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
 371         AliESDtrack iotrack;
 372         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 373         iotrack.SetTPCPoints(pt->GetPoints());
 374         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 375         //iotrack.SetTPCindex(i);
 376         fEvent->AddTrack(&iotrack);
 377         continue;
 378       }
 379
 380       if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.55) {
 381         AliESDtrack iotrack;
 382         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 383         iotrack.SetTPCPoints(pt->GetPoints());
 384         //iotrack.SetTPCindex(i);
 385         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 386         fEvent->AddTrack(&iotrack);
 387         continue;
 388       }
 389       //
 390       // short tracks  - maybe decays
 391
 392       if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.70) {
 393         Int_t found,foundable,shared;
 394         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
 395         if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2)){
 396           AliESDtrack iotrack;
 397           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 398           //iotrack.SetTPCindex(i);
 399           iotrack.SetTPCPoints(pt->GetPoints());
 400           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 401           fEvent->AddTrack(&iotrack);
 402           continue;
 403         }
 404       }
 405
 406       if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->fNFoundable))>0.8) {
 407         Int_t found,foundable,shared;
 408         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
 409         if (found<20) continue;
 410         if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
 411         //
 412         AliESDtrack iotrack;
 413         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 414         iotrack.SetTPCPoints(pt->GetPoints());
 415         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 416         //iotrack.SetTPCindex(i);
 417         fEvent->AddTrack(&iotrack);
 418         continue;
 419       }
 420       // short tracks  - secondaties
 421       //
 422       if ( (pt->GetNumberOfClusters()>30) ) {
 423         Int_t found,foundable,shared;
 424         pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
 425         if ( (found>20) && (pt->fNShared/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
 426           AliESDtrack iotrack;
 427           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 428           iotrack.SetTPCPoints(pt->GetPoints());
 429           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 430           //iotrack.SetTPCindex(i);
 431           fEvent->AddTrack(&iotrack);
 432           continue;
 433         }
 434       }
 435
 436       if ( (pt->GetNumberOfClusters()>15)) {
 437         Int_t found,foundable,shared;
 438         pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
 439         if (found<15) continue;
 440         if (pt->fNShared/float(pt->GetNumberOfClusters())>0.2) continue;
 441         if (float(found)/float(foundable)<0.8) continue;
 442         //
 443         AliESDtrack iotrack;
 444         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 445         iotrack.SetTPCPoints(pt->GetPoints());
 446         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 447         //iotrack.SetTPCindex(i);
 448         fEvent->AddTrack(&iotrack);
 449         continue;
 450       }
 451     }
 452   }
 453   printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
 454 }
 455
 456 void AliTPCtrackerMI::WriteTracks(TTree * tree)
 457 {
 458   //
 459   // write tracks from seed array to selected tree
 460   //
 461   fOutput  = tree;
 462   if (fOutput){
 463     AliTPCtrack *iotrack= new AliTPCtrack;
 464     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 465   }
 466   WriteTracks();
 467 }
 468
 469 void AliTPCtrackerMI::WriteTracks()
 470 {
 471   //
 472   // write tracks to the given output tree -
 473   // output specified with SetIO routine
 474   if (!fSeeds)  return;
 475   if (!fOutput){
 476     SetIO();
 477   }
 478
 479   if (fOutput){
 480     AliTPCtrack *iotrack= 0;
 481     Int_t nseed=fSeeds->GetEntriesFast();
 482     //for (Int_t i=0; i<nseed; i++) {
 483     //  iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
 484     //  if (iotrack) break;
 485     //}
 486     //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 487     TBranch * br = fOutput->GetBranch("tracks");
 488     br->SetAddress(&iotrack);
 489     //
 490     for (Int_t i=0; i<nseed; i++) {
 491       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
 492       if (!pt) continue;
 493       AliTPCtrack * track = new AliTPCtrack(*pt);
 494       iotrack = track;
 495       pt->fLab2 =i;
 496       //      br->SetAddress(&iotrack);
 497       fOutput->Fill();
 498       delete track;
 499       iotrack =0;
 500     }
 501     //fOutput->GetDirectory()->cd();
 502     //fOutput->Write();
 503   }
 504   // delete iotrack;
 505   //
 506   if (fSeedTree){
 507     //write the full seed information if specified in debug mode
 508
 509     AliTPCseed * vseed = new AliTPCseed;
 510     //
 511     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
 512     arrtr->ExpandCreateFast(160);
 513     //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
 514     //arrcl->ExpandCreateFast(160);
 515     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
 516     //
 517     vseed->fPoints = arrtr;
 518     vseed->fEPoints = arre;
 519     //    vseed->fClusterPoints = arrcl;
 520     //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
 521     TBranch * brseed = fSeedTree->GetBranch("seeds");
 522
 523     Int_t nseed=fSeeds->GetEntriesFast();
 524
 525     for (Int_t i=0; i<nseed; i++) {
 526       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
 527       if (!pt) continue;
 528       pt->fPoints = arrtr;
 529       //      pt->fClusterPoints = arrcl;
 530       pt->fEPoints       = arre;
 531       pt->RebuildSeed();
 532       vseed = pt;
 533       brseed->SetAddress(&vseed);
 534       fSeedTree->Fill();
 535       pt->fPoints  = 0;
 536       pt->fEPoints = 0;
 537       //      pt->fClusterPoints = 0;
 538     }
 539     fSeedTree->Write();
 540     if (fTreeDebug) fTreeDebug->Write();
 541   }
 542
 543 }
 544
 545
 546
 547
 548 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
 549   //
 550   //
 551   //seed->SetErrorY2(0.1);
 552   //return 0.1;
 553   //calculate look-up table at the beginning
 554   static Bool_t  ginit = kFALSE;
 555   static Float_t gnoise1,gnoise2,gnoise3;
 556   static Float_t ggg1[10000];
 557   static Float_t ggg2[10000];
 558   static Float_t ggg3[10000];
 559   static Float_t glandau1[10000];
 560   static Float_t glandau2[10000];
 561   static Float_t glandau3[10000];
 562   //
 563   static Float_t gcor01[500];
 564   static Float_t gcor02[500];
 565   static Float_t gcorp[500];
 566   //
 567
 568   //
 569   if (ginit==kFALSE){
 570     for (Int_t i=1;i<500;i++){
 571       Float_t rsigma = float(i)/100.;
 572       gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
 573       gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
 574       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 575     }
 576
 577     //
 578     for (Int_t i=3;i<10000;i++){
 579       //
 580       //
 581       // inner sector
 582       Float_t amp = float(i);
 583       Float_t padlength =0.75;
 584       gnoise1 = 0.0004/padlength;
 585       Float_t nel     = 0.268*amp;
 586       Float_t nprim   = 0.155*amp;
 587       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 588       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 589       if (glandau1[i]>1) glandau1[i]=1;
 590       glandau1[i]*=padlength*padlength/12.;
 591       //
 592       // outer short
 593       padlength =1.;
 594       gnoise2   = 0.0004/padlength;
 595       nel       = 0.3*amp;
 596       nprim     = 0.133*amp;
 597       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 598       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 599       if (glandau2[i]>1) glandau2[i]=1;
 600       glandau2[i]*=padlength*padlength/12.;
 601       //
 602       //
 603       // outer long
 604       padlength =1.5;
 605       gnoise3   = 0.0004/padlength;
 606       nel       = 0.3*amp;
 607       nprim     = 0.133*amp;
 608       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 609       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 610       if (glandau3[i]>1) glandau3[i]=1;
 611       glandau3[i]*=padlength*padlength/12.;
 612       //
 613     }
 614     ginit = kTRUE;
 615   }
 616   //
 617   //
 618   //
 619   Int_t amp = int(TMath::Abs(cl->GetQ()));
 620   if (amp>9999) {
 621     seed->SetErrorY2(1.);
 622     return 1.;
 623   }
 624   Float_t snoise2;
 625   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 626   Int_t ctype = cl->GetType();
 627   Float_t padlength= GetPadPitchLength(seed->fRow);
 628   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 629   angle2 = angle2/(1-angle2);
 630   //
 631   //cluster "quality"
 632   Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
 633   Float_t res;
 634   //
 635   if (fSectors==fInnerSec){
 636     snoise2 = gnoise1;
 637     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 638     if (ctype==0) res *= gcor01[rsigmay];
 639     if ((ctype>0)){
 640       res+=0.002;
 641       res*= gcorp[rsigmay];
 642     }
 643   }
 644   else {
 645     if (padlength<1.1){
 646       snoise2 = gnoise2;
 647       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 648       if (ctype==0) res *= gcor02[rsigmay];
 649       if ((ctype>0)){
 650         res+=0.002;
 651         res*= gcorp[rsigmay];
 652       }
 653     }
 654     else{
 655       snoise2 = gnoise3;
 656       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 657       if (ctype==0) res *= gcor02[rsigmay];
 658       if ((ctype>0)){
 659         res+=0.002;
 660         res*= gcorp[rsigmay];
 661       }
 662     }
 663   }
 664
 665   if (ctype<0){
 666     res+=0.005;
 667     res*=2.4;  // overestimate error 2 times
 668   }
 669   res+= snoise2;
 670
 671   if (res<2*snoise2)
 672     res = 2*snoise2;
 673
 674   seed->SetErrorY2(res);
 675   return res;
 676
 677
 678 }
 679
 680
 681
 682 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 683   //
 684   //
 685   //seed->SetErrorY2(0.1);
 686   //return 0.1;
 687   //calculate look-up table at the beginning
 688   static Bool_t  ginit = kFALSE;
 689   static Float_t gnoise1,gnoise2,gnoise3;
 690   static Float_t ggg1[10000];
 691   static Float_t ggg2[10000];
 692   static Float_t ggg3[10000];
 693   static Float_t glandau1[10000];
 694   static Float_t glandau2[10000];
 695   static Float_t glandau3[10000];
 696   //
 697   static Float_t gcor01[1000];
 698   static Float_t gcor02[1000];
 699   static Float_t gcorp[1000];
 700   //
 701
 702   //
 703   if (ginit==kFALSE){
 704     for (Int_t i=1;i<1000;i++){
 705       Float_t rsigma = float(i)/100.;
 706       gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
 707       gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
 708       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 709     }
 710
 711     //
 712     for (Int_t i=3;i<10000;i++){
 713       //
 714       //
 715       // inner sector
 716       Float_t amp = float(i);
 717       Float_t padlength =0.75;
 718       gnoise1 = 0.0004/padlength;
 719       Float_t nel     = 0.268*amp;
 720       Float_t nprim   = 0.155*amp;
 721       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 722       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 723       if (glandau1[i]>1) glandau1[i]=1;
 724       glandau1[i]*=padlength*padlength/12.;
 725       //
 726       // outer short
 727       padlength =1.;
 728       gnoise2   = 0.0004/padlength;
 729       nel       = 0.3*amp;
 730       nprim     = 0.133*amp;
 731       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 732       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 733       if (glandau2[i]>1) glandau2[i]=1;
 734       glandau2[i]*=padlength*padlength/12.;
 735       //
 736       //
 737       // outer long
 738       padlength =1.5;
 739       gnoise3   = 0.0004/padlength;
 740       nel       = 0.3*amp;
 741       nprim     = 0.133*amp;
 742       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 743       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 744       if (glandau3[i]>1) glandau3[i]=1;
 745       glandau3[i]*=padlength*padlength/12.;
 746       //
 747     }
 748     ginit = kTRUE;
 749   }
 750   //
 751   //
 752   //
 753   Int_t amp = int(TMath::Abs(cl->GetQ()));
 754   if (amp>9999) {
 755     seed->SetErrorY2(1.);
 756     return 1.;
 757   }
 758   Float_t snoise2;
 759   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 760   Int_t ctype = cl->GetType();
 761   Float_t padlength= GetPadPitchLength(seed->fRow);
 762   //
 763   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 764   //  if (angle2<0.6) angle2 = 0.6;
 765   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
 766   //
 767   //cluster "quality"
 768   Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
 769   Float_t res;
 770   //
 771   if (fSectors==fInnerSec){
 772     snoise2 = gnoise1;
 773     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 774     if (ctype==0) res *= gcor01[rsigmaz];
 775     if ((ctype>0)){
 776       res+=0.002;
 777       res*= gcorp[rsigmaz];
 778     }
 779   }
 780   else {
 781     if (padlength<1.1){
 782       snoise2 = gnoise2;
 783       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 784       if (ctype==0) res *= gcor02[rsigmaz];
 785       if ((ctype>0)){
 786         res+=0.002;
 787         res*= gcorp[rsigmaz];
 788       }
 789     }
 790     else{
 791       snoise2 = gnoise3;
 792       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 793       if (ctype==0) res *= gcor02[rsigmaz];
 794       if ((ctype>0)){
 795         res+=0.002;
 796         res*= gcorp[rsigmaz];
 797       }
 798     }
 799   }
 800
 801   if (ctype<0){
 802     res+=0.002;
 803     res*=1.3;
 804   }
 805   if ((ctype<0) &&amp<70){
 806     res+=0.002;
 807     res*=1.3;
 808   }
 809   res += snoise2;
 810   if (res<2*snoise2)
 811      res = 2*snoise2;
 812   if (res>3) res =3;
 813   seed->SetErrorZ2(res);
 814   return res;
 815 }
 816
 817
 818
 819 /*
 820 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 821   //
 822   //
 823   //seed->SetErrorZ2(0.1);
 824   //return 0.1;
 825
 826   Float_t snoise2;
 827   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 828   //
 829   Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
 830   Int_t ctype = cl->GetType();
 831   Float_t amp = TMath::Abs(cl->GetQ());
 832
 833   Float_t nel;
 834   Float_t nprim;
 835   //
 836   Float_t landau=2 ;    //landau fluctuation part
 837   Float_t gg=2;         // gg fluctuation part
 838   Float_t padlength= GetPadPitchLength(seed->GetX());
 839
 840   if (fSectors==fInnerSec){
 841     snoise2 = 0.0004/padlength;
 842     nel     = 0.268*amp;
 843     nprim   = 0.155*amp;
 844     gg      = (2+0.001*nel/(padlength*padlength))/nel;
 845     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 846     if (landau>1) landau=1;
 847   }
 848   else {
 849     snoise2 = 0.0004/padlength;
 850     nel     = 0.3*amp;
 851     nprim   = 0.133*amp;
 852     gg      = (2+0.0008*nel/(padlength*padlength))/nel;
 853     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 854     if (landau>1) landau=1;
 855   }
 856   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
 857
 858   //
 859   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 860   angle2 = TMath::Sqrt((1-angle2));
 861   if (angle2<0.6) angle2 = 0.6;
 862   //angle2 = 1;
 863
 864   Float_t angle = seed->GetTgl()/angle2;
 865   Float_t angular = landau*angle*angle*padlength*padlength/12.;
 866   Float_t res = sdiff + angular;
 867
 868
 869   if ((ctype==0) && (fSectors ==fOuterSec))
 870     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
 871
 872   if ((ctype==0) && (fSectors ==fInnerSec))
 873     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
 874
 875   if ((ctype>0)){
 876     res+=0.005;
 877     res*= TMath::Power(rsigmaz+0.5,1.5);  //0.31+0.147*ctype;
 878   }
 879   if (ctype<0){
 880     res+=0.002;
 881     res*=1.3;
 882   }
 883   if ((ctype<0) &&amp<70){
 884     res+=0.002;
 885     res*=1.3;
 886   }
 887   res += snoise2;
 888   if (res<2*snoise2)
 889      res = 2*snoise2;
 890
 891   seed->SetErrorZ2(res);
 892   return res;
 893 }
 894 */
 895
 896
 897
 898 void AliTPCseed::Reset(Bool_t all)
 899 {
 900   //
 901   //
 902   SetNumberOfClusters(0);
 903   fNFoundable = 0;
 904   SetChi2(0);
 905   ResetCovariance();
 906   /*
 907   if (fTrackPoints){
 908     for (Int_t i=0;i<8;i++){
 909       delete [] fTrackPoints[i];
 910     }
 911     delete fTrackPoints;
 912     fTrackPoints =0;
 913   }
 914   */
 915
 916   if (all){
 917     for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
 918     for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 919   }
 920
 921 }
 922
 923
 924 void AliTPCseed::Modify(Double_t factor)
 925 {
 926
 927   //------------------------------------------------------------------
 928   //This function makes a track forget its history :)
 929   //------------------------------------------------------------------
 930   if (factor<=0) {
 931     ResetCovariance();
 932     return;
 933   }
 934   fC00*=factor;
 935   fC10*=0;  fC11*=factor;
 936   fC20*=0;  fC21*=0;  fC22*=factor;
 937   fC30*=0;  fC31*=0;  fC32*=0;  fC33*=factor;
 938   fC40*=0;  fC41*=0;  fC42*=0;  fC43*=0;  fC44*=factor;
 939   SetNumberOfClusters(0);
 940   fNFoundable =0;
 941   SetChi2(0);
 942   fRemoval = 0;
 943   fCurrentSigmaY2 = 0.000005;
 944   fCurrentSigmaZ2 = 0.000005;
 945   fNoCluster     = 0;
 946   //fFirstPoint = 160;
 947   //fLastPoint  = 0;
 948 }
 949
 950
 951
 952
 953 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
 954 {
 955   //-----------------------------------------------------------------
 956   // This function find proloncation of a track to a reference plane x=xk.
 957   // doesn't change internal state of the track
 958   //-----------------------------------------------------------------
 959
 960   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
 961
 962   if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
 963     return 0;
 964   }
 965
 966   //  Double_t y1=fP0, z1=fP1;
 967   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
 968   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
 969
 970   y = fP0;
 971   z = fP1;
 972   //y += dx*(c1+c2)/(r1+r2);
 973   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 974
 975   Double_t dy = dx*(c1+c2)/(r1+r2);
 976   Double_t dz = 0;
 977   //
 978   Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
 979   /*
 980   if (TMath::Abs(delta)>0.0001){
 981     dz = fP3*TMath::ASin(delta)/fP4;
 982   }else{
 983     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
 984   }
 985   */
 986   dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
 987   //
 988   y+=dy;
 989   z+=dz;
 990
 991
 992   return 1;
 993 }
 994
 995
 996 //_____________________________________________________________________________
 997 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
 998 {
 999   //-----------------------------------------------------------------
1000   // This function calculates a predicted chi2 increment.
1001   //-----------------------------------------------------------------
1002   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
1003   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1004   r00+=fC00; r01+=fC10; r11+=fC11;
1005
1006   Double_t det=r00*r11 - r01*r01;
1007   if (TMath::Abs(det) < 1.e-10) {
1008     Int_t n=GetNumberOfClusters();
1009     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
1010     return 1e10;
1011   }
1012   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1013
1014   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1015
1016   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
1017 }
1018
1019
1020 //_________________________________________________________________________________________
1021
1022
1023 Int_t AliTPCseed::Compare(const TObject *o) const {
1024   //-----------------------------------------------------------------
1025   // This function compares tracks according to the sector - for given sector according z
1026   //-----------------------------------------------------------------
1027   AliTPCseed *t=(AliTPCseed*)o;
1028
1029   if (fSort == 0){
1030     if (t->fRelativeSector>fRelativeSector) return -1;
1031     if (t->fRelativeSector<fRelativeSector) return 1;
1032     Double_t z2 = t->GetZ();
1033     Double_t z1 = GetZ();
1034     if (z2>z1) return 1;
1035     if (z2<z1) return -1;
1036     return 0;
1037   }
1038   else {
1039     Float_t f2 =1;
1040     f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
1041     if (t->fBConstrain) f2=1.2;
1042
1043     Float_t f1 =1;
1044     f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
1045
1046     if (fBConstrain)   f1=1.2;
1047
1048     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
1049     else return +1;
1050   }
1051 }
1052
1053 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
1054 {
1055   //rotate to track "local coordinata
1056   Float_t x = seed->GetX();
1057   Float_t y = seed->GetY();
1058   Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
1059
1060   if (y > ymax) {
1061     seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
1062     if (!seed->Rotate(fSectors->GetAlpha()))
1063       return;
1064   } else if (y <-ymax) {
1065     seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
1066     if (!seed->Rotate(-fSectors->GetAlpha()))
1067       return;
1068   }
1069
1070 }
1071
1072
1073
1074
1075 //_____________________________________________________________________________
1076 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
1077   //-----------------------------------------------------------------
1078   // This function associates a cluster with this track.
1079   //-----------------------------------------------------------------
1080   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
1081
1082   r00+=fC00; r01+=fC10; r11+=fC11;
1083   Double_t det=r00*r11 - r01*r01;
1084   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
1085
1086   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
1087   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
1088   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
1089   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
1090   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
1091
1092   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
1093   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
1094   if (TMath::Abs(cur*fX-eta) >= 0.9) {
1095     return 0;
1096   }
1097
1098   fP0 += k00*dy + k01*dz;
1099   fP1 += k10*dy + k11*dz;
1100   fP2  = eta;
1101   fP3 += k30*dy + k31*dz;
1102   fP4  = cur;
1103
1104   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1105   Double_t c12=fC21, c13=fC31, c14=fC41;
1106
1107   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1108   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
1109   fC40-=k00*c04+k01*c14;
1110
1111   fC11-=k10*c01+k11*fC11;
1112   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
1113   fC41-=k10*c04+k11*c14;
1114
1115   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
1116   fC42-=k20*c04+k21*c14;
1117
1118   fC33-=k30*c03+k31*c13;
1119   fC43-=k40*c03+k41*c13;
1120
1121   fC44-=k40*c04+k41*c14;
1122
1123   Int_t n=GetNumberOfClusters();
1124   //  fIndex[n]=index;
1125   SetNumberOfClusters(n+1);
1126   SetChi2(GetChi2()+chisq);
1127
1128   return 1;
1129 }
1130
1131
1132
1133 //_____________________________________________________________________________
1134 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
1135                    Double_t x2,Double_t y2,
1136                    Double_t x3,Double_t y3)
1137 {
1138   //-----------------------------------------------------------------
1139   // Initial approximation of the track curvature
1140   //-----------------------------------------------------------------
1141   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1142   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1143                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1144   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1145                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1146
1147   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1148   if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1149   return -xr*yr/sqrt(xr*xr+yr*yr);
1150 }
1151
1152
1153
1154 //_____________________________________________________________________________
1155 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
1156                    Double_t x2,Double_t y2,
1157                    Double_t x3,Double_t y3)
1158 {
1159   //-----------------------------------------------------------------
1160   // Initial approximation of the track curvature
1161   //-----------------------------------------------------------------
1162   x3 -=x1;
1163   x2 -=x1;
1164   y3 -=y1;
1165   y2 -=y1;
1166   //
1167   Double_t det = x3*y2-x2*y3;
1168   if (det==0) {
1169     return 100;
1170   }
1171   //
1172   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1173   Double_t x0 = x3*0.5-y3*u;
1174   Double_t y0 = y3*0.5+x3*u;
1175   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1176   if (det<0) c2*=-1;
1177   return c2;
1178 }
1179
1180
1181 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
1182                    Double_t x2,Double_t y2,
1183                    Double_t x3,Double_t y3)
1184 {
1185   //-----------------------------------------------------------------
1186   // Initial approximation of the track curvature
1187   //-----------------------------------------------------------------
1188   x3 -=x1;
1189   x2 -=x1;
1190   y3 -=y1;
1191   y2 -=y1;
1192   //
1193   Double_t det = x3*y2-x2*y3;
1194   if (det==0) {
1195     return 100;
1196   }
1197   //
1198   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1199   Double_t x0 = x3*0.5-y3*u;
1200   Double_t y0 = y3*0.5+x3*u;
1201   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1202   if (det<0) c2*=-1;
1203   x0+=x1;
1204   x0*=c2;
1205   return x0;
1206 }
1207
1208
1209
1210 //_____________________________________________________________________________
1211 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1212                    Double_t x2,Double_t y2,
1213                    Double_t x3,Double_t y3)
1214 {
1215   //-----------------------------------------------------------------
1216   // Initial approximation of the track curvature times center of curvature
1217   //-----------------------------------------------------------------
1218   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1219   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1220                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1221   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1222                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1223
1224   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1225
1226   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1227 }
1228
1229 //_____________________________________________________________________________
1230 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1231                    Double_t x2,Double_t y2,
1232                    Double_t z1,Double_t z2)
1233 {
1234   //-----------------------------------------------------------------
1235   // Initial approximation of the tangent of the track dip angle
1236   //-----------------------------------------------------------------
1237   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1238 }
1239
1240
1241 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1242                    Double_t x2,Double_t y2,
1243                    Double_t z1,Double_t z2, Double_t c)
1244 {
1245   //-----------------------------------------------------------------
1246   // Initial approximation of the tangent of the track dip angle
1247   //-----------------------------------------------------------------
1248
1249   //  Double_t angle1;
1250
1251   //angle1    =  (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1252   //
1253   Double_t d  =  TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1254   if (TMath::Abs(d*c*0.5)>1) return 0;
1255   //  Double_t   angle2    =  TMath::ASin(d*c*0.5);
1256   //  Double_t   angle2    =  AliTPCFastMath::FastAsin(d*c*0.5);
1257   Double_t   angle2    = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1258
1259   angle2  = (z1-z2)*c/(angle2*2.);
1260   return angle2;
1261 }
1262
1263 Bool_t   AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1264 {//-----------------------------------------------------------------
1265   // This function find proloncation of a track to a reference plane x=x2.
1266   //-----------------------------------------------------------------
1267
1268   Double_t dx=x2-x1;
1269
1270   if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1271     return kFALSE;
1272   }
1273
1274   Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1275   Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1276   y = x[0];
1277   z = x[1];
1278
1279   Double_t dy = dx*(c1+c2)/(r1+r2);
1280   Double_t dz = 0;
1281   //
1282   Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1283
1284   if (TMath::Abs(delta)>0.01){
1285     dz = x[3]*TMath::ASin(delta)/x[4];
1286   }else{
1287     dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1288   }
1289
1290   //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1291
1292   y+=dy;
1293   z+=dz;
1294
1295   return kTRUE;
1296 }
1297
1298 Int_t  AliTPCtrackerMI::LoadClusters (TTree *tree)
1299 {
1300   //
1301   //
1302   fInput = tree;
1303   return LoadClusters();
1304 }
1305
1306 Int_t  AliTPCtrackerMI::LoadClusters()
1307 {
1308   //
1309   // load clusters to the memory
1310   AliTPCClustersRow *clrow= new AliTPCClustersRow;
1311   clrow->SetClass("AliTPCclusterMI");
1312   clrow->SetArray(0);
1313   clrow->GetArray()->ExpandCreateFast(10000);
1314   //
1315   //  TTree * tree = fClustersArray.GetTree();
1316
1317   TTree * tree = fInput;
1318   TBranch * br = tree->GetBranch("Segment");
1319   br->SetAddress(&clrow);
1320   //
1321   Int_t j=Int_t(tree->GetEntries());
1322   for (Int_t i=0; i<j; i++) {
1323     br->GetEntry(i);
1324     //
1325     Int_t sec,row;
1326     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1327     //
1328     AliTPCRow * tpcrow=0;
1329     Int_t left=0;
1330     if (sec<fkNIS*2){
1331       tpcrow = &(fInnerSec[sec%fkNIS][row]);
1332       left = sec/fkNIS;
1333     }
1334     else{
1335       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1336       left = (sec-fkNIS*2)/fkNOS;
1337     }
1338     if (left ==0){
1339       tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1340       tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1341       for (Int_t i=0;i<tpcrow->fN1;i++)
1342         tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1343     }
1344     if (left ==1){
1345       tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1346       tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1347       for (Int_t i=0;i<tpcrow->fN2;i++)
1348         tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1349     }
1350   }
1351   //
1352   delete clrow;
1353   LoadOuterSectors();
1354   LoadInnerSectors();
1355   return 0;
1356 }
1357
1358
1359 void AliTPCtrackerMI::UnloadClusters()
1360 {
1361   //
1362   // unload clusters from the memory
1363   //
1364   Int_t nrows = fOuterSec->GetNRows();
1365   for (Int_t sec = 0;sec<fkNOS;sec++)
1366     for (Int_t row = 0;row<nrows;row++){
1367       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1368       //      if (tpcrow){
1369       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1370       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1371       //}
1372       tpcrow->ResetClusters();
1373     }
1374   //
1375   nrows = fInnerSec->GetNRows();
1376   for (Int_t sec = 0;sec<fkNIS;sec++)
1377     for (Int_t row = 0;row<nrows;row++){
1378       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1379       //if (tpcrow){
1380       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1381       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1382       //}
1383       tpcrow->ResetClusters();
1384     }
1385
1386   return ;
1387 }
1388
1389
1390 //_____________________________________________________________________________
1391 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1392   //-----------------------------------------------------------------
1393   // This function fills outer TPC sectors with clusters.
1394   //-----------------------------------------------------------------
1395   Int_t nrows = fOuterSec->GetNRows();
1396   UInt_t index=0;
1397   for (Int_t sec = 0;sec<fkNOS;sec++)
1398     for (Int_t row = 0;row<nrows;row++){
1399       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1400       Int_t sec2 = sec+2*fkNIS;
1401       //left
1402       Int_t ncl = tpcrow->fN1;
1403       while (ncl--) {
1404         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1405         index=(((sec2<<8)+row)<<16)+ncl;
1406         tpcrow->InsertCluster(c,index);
1407       }
1408       //right
1409       ncl = tpcrow->fN2;
1410       while (ncl--) {
1411         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1412         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1413         tpcrow->InsertCluster(c,index);
1414       }
1415       //
1416       // write indexes for fast acces
1417       //
1418       for (Int_t i=0;i<510;i++)
1419         tpcrow->fFastCluster[i]=-1;
1420       for (Int_t i=0;i<tpcrow->GetN();i++){
1421         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1422         tpcrow->fFastCluster[zi]=i;  // write index
1423       }
1424       Int_t last = 0;
1425       for (Int_t i=0;i<510;i++){
1426         if (tpcrow->fFastCluster[i]<0)
1427           tpcrow->fFastCluster[i] = last;
1428         else
1429           last = tpcrow->fFastCluster[i];
1430       }
1431     }
1432   fN=fkNOS;
1433   fSectors=fOuterSec;
1434   return 0;
1435 }
1436
1437
1438 //_____________________________________________________________________________
1439 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1440   //-----------------------------------------------------------------
1441   // This function fills inner TPC sectors with clusters.
1442   //-----------------------------------------------------------------
1443   Int_t nrows = fInnerSec->GetNRows();
1444   UInt_t index=0;
1445   for (Int_t sec = 0;sec<fkNIS;sec++)
1446     for (Int_t row = 0;row<nrows;row++){
1447       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1448       //
1449       //left
1450       Int_t ncl = tpcrow->fN1;
1451       while (ncl--) {
1452         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1453         index=(((sec<<8)+row)<<16)+ncl;
1454         tpcrow->InsertCluster(c,index);
1455       }
1456       //right
1457       ncl = tpcrow->fN2;
1458       while (ncl--) {
1459         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1460         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1461         tpcrow->InsertCluster(c,index);
1462       }
1463       //
1464       // write indexes for fast acces
1465       //
1466       for (Int_t i=0;i<510;i++)
1467         tpcrow->fFastCluster[i]=-1;
1468       for (Int_t i=0;i<tpcrow->GetN();i++){
1469         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1470         tpcrow->fFastCluster[zi]=i;  // write index
1471       }
1472       Int_t last = 0;
1473       for (Int_t i=0;i<510;i++){
1474         if (tpcrow->fFastCluster[i]<0)
1475           tpcrow->fFastCluster[i] = last;
1476         else
1477           last = tpcrow->fFastCluster[i];
1478       }
1479
1480     }
1481
1482   fN=fkNIS;
1483   fSectors=fInnerSec;
1484   return 0;
1485 }
1486
1487
1488
1489 //_________________________________________________________________________
1490 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1491   //--------------------------------------------------------------------
1492   //       Return pointer to a given cluster
1493   //--------------------------------------------------------------------
1494   Int_t sec=(index&0xff000000)>>24;
1495   Int_t row=(index&0x00ff0000)>>16;
1496   Int_t ncl=(index&0x00007fff)>>00;
1497
1498   const AliTPCRow * tpcrow=0;
1499   AliTPCclusterMI * clrow =0;
1500
1501   if (sec<fkNIS*2){
1502     tpcrow = &(fInnerSec[sec%fkNIS][row]);
1503     if (tpcrow==0) return 0;
1504
1505     if (sec<fkNIS) {
1506       if (tpcrow->fN1<=ncl) return 0;
1507       clrow = tpcrow->fClusters1;
1508     }
1509     else {
1510       if (tpcrow->fN2<=ncl) return 0;
1511       clrow = tpcrow->fClusters2;
1512     }
1513   }
1514   else {
1515     tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1516     if (tpcrow==0) return 0;
1517
1518     if (sec-2*fkNIS<fkNOS) {
1519       if (tpcrow->fN1<=ncl) return 0;
1520       clrow = tpcrow->fClusters1;
1521     }
1522     else {
1523       if (tpcrow->fN2<=ncl) return 0;
1524       clrow = tpcrow->fClusters2;
1525     }
1526   }
1527
1528   return &(clrow[ncl]);
1529
1530 }
1531
1532
1533
1534 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1535   //-----------------------------------------------------------------
1536   // This function tries to find a track prolongation to next pad row
1537   //-----------------------------------------------------------------
1538   //
1539   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1540   AliTPCclusterMI *cl=0;
1541   Int_t tpcindex= t.GetClusterIndex2(nr);
1542   //
1543   // update current shape info every 5 pad-row
1544   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1545     GetShape(&t,nr);
1546     //}
1547   //
1548   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters
1549     //
1550     if (tpcindex==-1) return 0; //track in dead zone
1551     if (tpcindex>0){     //
1552       cl = t.fClusterPointer[nr];
1553       if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1554       t.fCurrentClusterIndex1 = tpcindex;
1555     }
1556     if (cl){
1557       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
1558       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1559       //
1560       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1561       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1562
1563       if (TMath::Abs(angle-t.GetAlpha())>0.001){
1564         Double_t rotation = angle-t.GetAlpha();
1565         t.fRelativeSector= relativesector;
1566         t.Rotate(rotation);
1567       }
1568       t.PropagateTo(x);
1569       //
1570       t.fCurrentCluster = cl;
1571       t.fRow = nr;
1572       Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1573       if ((tpcindex&0x8000)==0) accept =0;
1574       if (accept<3) {
1575         //if founded cluster is acceptible
1576         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
1577           t.fErrorY2 += 0.03;
1578           t.fErrorZ2 += 0.03;
1579           t.fErrorY2 *= 3;
1580           t.fErrorZ2 *= 3;
1581         }
1582         t.fNFoundable++;
1583         UpdateTrack(&t,accept);
1584         return 1;
1585       }
1586     }
1587   }
1588   if (fIteration>1) return 0;  // not look for new cluster during refitting
1589   //
1590   UInt_t index=0;
1591   if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;
1592   Double_t  y=t.GetYat(x);
1593   if (TMath::Abs(y)>ymax){
1594     if (y > ymax) {
1595       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1596       if (!t.Rotate(fSectors->GetAlpha()))
1597         return 0;
1598     } else if (y <-ymax) {
1599       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1600       if (!t.Rotate(-fSectors->GetAlpha()))
1601         return 0;
1602     }
1603     //return 1;
1604   }
1605   //
1606   if (!t.PropagateTo(x)) {
1607     if (fIteration==0) t.fRemoval = 10;
1608     return 0;
1609   }
1610   y=t.GetY();
1611   Double_t z=t.GetZ();
1612   //
1613   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1614   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1615   Double_t  roady  =1.;
1616   Double_t  roadz = 1.;
1617   //
1618   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1619     t.fInDead = kTRUE;
1620     t.SetClusterIndex2(nr,-1);
1621     return 0;
1622   }
1623   else
1624     {
1625       if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10)) t.fNFoundable++;
1626       else
1627         return 0;
1628     }
1629   //calculate
1630   if (krow) {
1631     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1632     cl = krow.FindNearest2(y,z,roady,roadz,index);
1633     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1634   }
1635   if (cl) {
1636     t.fCurrentCluster = cl;
1637     t.fRow = nr;
1638     if (fIteration==2&&cl->IsUsed(10)) return 0;
1639     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1640     if (fIteration==2&&cl->IsUsed(11)) {
1641       t.fErrorY2 += 0.03;
1642       t.fErrorZ2 += 0.03;
1643       t.fErrorY2 *= 3;
1644       t.fErrorZ2 *= 3;
1645     }
1646     /*
1647     if (t.fCurrentCluster->IsUsed(10)){
1648       //
1649       //
1650
1651       t.fNShared++;
1652       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1653         t.fRemoval =10;
1654         return 0;
1655       }
1656     }
1657     */
1658     if (accept<3) UpdateTrack(&t,accept);
1659
1660   } else {
1661     if ( fIteration==0 && t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1662
1663   }
1664   return 1;
1665 }
1666
1667 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1668   //-----------------------------------------------------------------
1669   // This function tries to find a track prolongation to next pad row
1670   //-----------------------------------------------------------------
1671   //
1672   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1673   Double_t y,z;
1674   if (!t.GetProlongation(x,y,z)) {
1675     t.fRemoval = 10;
1676     return 0;
1677   }
1678   //
1679   //
1680   if (TMath::Abs(y)>ymax){
1681
1682     if (y > ymax) {
1683       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1684       if (!t.Rotate(fSectors->GetAlpha()))
1685         return 0;
1686     } else if (y <-ymax) {
1687       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1688       if (!t.Rotate(-fSectors->GetAlpha()))
1689         return 0;
1690     }
1691     if (!t.PropagateTo(x)) {
1692       return 0;
1693     }
1694     t.GetProlongation(x,y,z);
1695   }
1696   //
1697   // update current shape info every 3 pad-row
1698   if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1699     //    t.fCurrentSigmaY = GetSigmaY(&t);
1700     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1701     GetShape(&t,nr);
1702   }
1703   //
1704   AliTPCclusterMI *cl=0;
1705   UInt_t index=0;
1706
1707
1708   //Int_t nr2 = nr;
1709   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1710   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1711   Double_t  roady  =1.;
1712   Double_t  roadz = 1.;
1713   //
1714   Int_t row = nr;
1715   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1716     t.fInDead = kTRUE;
1717     t.SetClusterIndex2(row,-1);
1718     return 0;
1719   }
1720   else
1721     {
1722       if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1723     }
1724   //calculate
1725
1726   if ((cl==0)&&(krow)) {
1727     //    cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1728     cl = krow.FindNearest2(y,z,roady,roadz,index);
1729
1730     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1731   }
1732
1733   if (cl) {
1734     t.fCurrentCluster = cl;
1735     //    Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1736     //if (accept<3){
1737       t.SetClusterIndex2(row,index);
1738       t.fClusterPointer[row] = cl;
1739       //}
1740   }
1741   return 1;
1742 }
1743
1744
1745
1746 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
1747   //-----------------------------------------------------------------
1748   // This function tries to find a track prolongation to next pad row
1749   //-----------------------------------------------------------------
1750   t.fCurrentCluster  = 0;
1751   t.fCurrentClusterIndex1 = 0;
1752
1753   Double_t xt=t.GetX();
1754   Int_t     row = GetRowNumber(xt)-1;
1755   Double_t  ymax= GetMaxY(nr);
1756
1757   if (row < nr) return 1; // don't prolongate if not information until now -
1758   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1759     t.fRemoval =10;
1760     return 0;  // not prolongate strongly inclined tracks
1761   }
1762   if (TMath::Abs(t.GetSnp())>0.95) {
1763     t.fRemoval =10;
1764     return 0;  // not prolongate strongly inclined tracks
1765   }
1766
1767   Double_t x= GetXrow(nr);
1768   Double_t y,z;
1769   //t.PropagateTo(x+0.02);
1770   //t.PropagateTo(x+0.01);
1771   if (!t.PropagateTo(x)){
1772     return 0;
1773   }
1774   //
1775   y=t.GetY();
1776   z=t.GetZ();
1777
1778   if (TMath::Abs(y)>ymax){
1779     if (y > ymax) {
1780       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1781       if (!t.Rotate(fSectors->GetAlpha()))
1782         return 0;
1783     } else if (y <-ymax) {
1784       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1785       if (!t.Rotate(-fSectors->GetAlpha()))
1786         return 0;
1787     }
1788     //    if (!t.PropagateTo(x)){
1789     //  return 0;
1790     //}
1791     return 1;
1792     //y = t.GetY();
1793   }
1794   //
1795
1796   AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1797
1798   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1799     t.fInDead = kTRUE;
1800     t.SetClusterIndex2(nr,-1);
1801     return 0;
1802   }
1803   else
1804     {
1805       if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10)) t.fNFoundable++;
1806       else
1807         return 0;
1808     }
1809
1810   // update current
1811   if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1812     //    t.fCurrentSigmaY = GetSigmaY(&t);
1813     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1814     GetShape(&t,nr);
1815   }
1816
1817   AliTPCclusterMI *cl=0;
1818   UInt_t index=0;
1819   //
1820   Double_t roady = 1.;
1821   Double_t roadz = 1.;
1822   //
1823
1824   if (!cl){
1825     index = t.GetClusterIndex2(nr);
1826     if ( (index>0) && (index&0x8000)==0){
1827       cl = t.fClusterPointer[nr];
1828       if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1829       t.fCurrentClusterIndex1 = index;
1830       if (cl) {
1831         t.fCurrentCluster  = cl;
1832         return 1;
1833       }
1834     }
1835   }
1836
1837   if (krow) {
1838     //cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1839     cl = krow.FindNearest2(y,z,roady,roadz,index);
1840   }
1841
1842   if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1843   t.fCurrentCluster  = cl;
1844
1845   return 1;
1846 }
1847
1848
1849 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1850   //-----------------------------------------------------------------
1851   // This function tries to find a track prolongation to next pad row
1852   //-----------------------------------------------------------------
1853
1854   //update error according neighborhoud
1855
1856   if (t.fCurrentCluster) {
1857     t.fRow = nr;
1858     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1859
1860     if (t.fCurrentCluster->IsUsed(10)){
1861       //
1862       //
1863       //  t.fErrorZ2*=2;
1864       //  t.fErrorY2*=2;
1865       t.fNShared++;
1866       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1867         t.fRemoval =10;
1868         return 0;
1869       }
1870     }
1871     if (fIteration>0) accept = 0;
1872     if (accept<3)  UpdateTrack(&t,accept);
1873
1874   } else {
1875     if (fIteration==0){
1876       if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1877       if (  t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1878
1879       if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1880     }
1881   }
1882   return 1;
1883 }
1884
1885
1886
1887 //_____________________________________________________________________________
1888 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1889   //-----------------------------------------------------------------
1890   // This function tries to find a track prolongation.
1891   //-----------------------------------------------------------------
1892   Double_t xt=t.GetX();
1893   //
1894   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1895   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1896   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1897   //
1898   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1899
1900   Int_t first = GetRowNumber(xt)-1;
1901   for (Int_t nr= first; nr>=rf; nr-=step) {
1902     // update kink info
1903     if (t.GetKinkIndexes()[0]>0){
1904       for (Int_t i=0;i<3;i++){
1905         Int_t index = t.GetKinkIndexes()[i];
1906         if (index==0) break;
1907         if (index<0) continue;
1908         //
1909         AliESDkink * kink = fEvent->GetKink(index-1);
1910         if (!kink){
1911           printf("PROBLEM\n");
1912         }
1913         else{
1914           Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1915           if (kinkrow==nr){
1916             AliExternalTrackParam paramd(t);
1917             kink->SetDaughter(paramd);
1918             kink->SetStatus(2,5);
1919             kink->Update();
1920           }
1921         }
1922       }
1923     }
1924
1925     if (nr==80) t.UpdateReference();
1926     if (nr<fInnerSec->GetNRows())
1927       fSectors = fInnerSec;
1928     else
1929       fSectors = fOuterSec;
1930     if (FollowToNext(t,nr)==0)
1931       if (!t.IsActive())
1932         return 0;
1933
1934   }
1935   return 1;
1936 }
1937
1938
1939 //_____________________________________________________________________________
1940 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1941   //-----------------------------------------------------------------
1942   // This function tries to find a track prolongation.
1943   //-----------------------------------------------------------------
1944   Double_t xt=t.GetX();
1945   //
1946   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1947   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1948   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1949   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1950
1951   for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1952
1953     if (FollowToNextFast(t,nr)==0)
1954       if (!t.IsActive()) return 0;
1955
1956   }
1957   return 1;
1958 }
1959
1960
1961
1962
1963
1964 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1965   //-----------------------------------------------------------------
1966   // This function tries to find a track prolongation.
1967   //-----------------------------------------------------------------
1968   //
1969   Double_t xt=t.GetX();
1970   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1971   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1972   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1973   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1974
1975   Int_t first = t.fFirstPoint;
1976   if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
1977   //
1978   if (first<0) first=0;
1979   for (Int_t nr=first; nr<=rf; nr++) {
1980     //if ( (t.GetSnp()<0.9))
1981     if (t.GetKinkIndexes()[0]<0){
1982       for (Int_t i=0;i<3;i++){
1983         Int_t index = t.GetKinkIndexes()[i];
1984         if (index==0) break;
1985         if (index>0) continue;
1986         index = TMath::Abs(index);
1987         AliESDkink * kink = fEvent->GetKink(index-1);
1988         if (!kink){
1989           printf("PROBLEM\n");
1990         }
1991         else{
1992           Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
1993           if (kinkrow==nr){
1994             AliExternalTrackParam paramm(t);
1995             kink->SetMother(paramm);
1996             kink->SetStatus(2,1);
1997             kink->Update();
1998           }
1999         }
2000       }
2001     }
2002     //
2003     if (nr<fInnerSec->GetNRows())
2004       fSectors = fInnerSec;
2005     else
2006       fSectors = fOuterSec;
2007     FollowToNext(t,nr);
2008   }
2009   return 1;
2010 }
2011
2012
2013
2014
2015
2016 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2017 {
2018   //
2019   //
2020   sum1=0;
2021   sum2=0;
2022   Int_t sum=0;
2023   //
2024   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2025   dz2*=dz2;
2026
2027   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2028   dy2*=dy2;
2029   Float_t distance = TMath::Sqrt(dz2+dy2);
2030   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2031
2032   //  Int_t offset =0;
2033   Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
2034   Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
2035   if (lastpoint>160)
2036     lastpoint =160;
2037   if (firstpoint<0)
2038     firstpoint = 0;
2039   if (firstpoint>lastpoint) {
2040     firstpoint =lastpoint;
2041     //    lastpoint  =160;
2042   }
2043
2044
2045   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2046     if (s1->GetClusterIndex2(i)>0) sum1++;
2047     if (s2->GetClusterIndex2(i)>0) sum2++;
2048     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2049       sum++;
2050     }
2051   }
2052   if (sum<5) return 0;
2053
2054   Float_t summin = TMath::Min(sum1+1,sum2+1);
2055   Float_t ratio = (sum+1)/Float_t(summin);
2056   return ratio;
2057 }
2058
2059 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2060 {
2061   //
2062   //
2063   if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
2064   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
2065
2066   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2067   dz2*=dz2;
2068   Float_t dy2 =(s1->GetY() - s2->GetY());
2069   dy2*=dy2;
2070   Float_t distance = dz2+dy2;
2071   if (distance>325.) return ; // if there are far away  - not overlap - to reduce combinatorics
2072
2073   //
2074   Int_t sumshared=0;
2075   //
2076   Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
2077   Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
2078   //
2079   if (firstpoint>=lastpoint-5) return;;
2080
2081   for (Int_t i=firstpoint;i<lastpoint;i++){
2082     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2083     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2084       sumshared++;
2085     }
2086   }
2087   if (sumshared>4){
2088     // sign clusters
2089     //
2090     for (Int_t i=firstpoint;i<lastpoint;i++){
2091       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2092       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2093         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2094         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);;
2095         if (s1->IsActive()&&s2->IsActive()){
2096           p1->fIsShared = kTRUE;
2097           p2->fIsShared = kTRUE;
2098         }
2099       }
2100     }
2101   }
2102   //
2103   if (sumshared>10){
2104     for (Int_t i=0;i<4;i++){
2105       if (s1->fOverlapLabels[3*i]==0){
2106         s1->fOverlapLabels[3*i] = s2->GetLabel();
2107         s1->fOverlapLabels[3*i+1] = sumshared;
2108         s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
2109         break;
2110       }
2111     }
2112     for (Int_t i=0;i<4;i++){
2113       if (s2->fOverlapLabels[3*i]==0){
2114         s2->fOverlapLabels[3*i] = s1->GetLabel();
2115         s2->fOverlapLabels[3*i+1] = sumshared;
2116         s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
2117         break;
2118       }
2119     }
2120   }
2121
2122 }
2123
2124 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
2125 {
2126   //
2127   //sort trackss according sectors
2128   //
2129   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2130     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2131     if (!pt) continue;
2132     //if (pt) RotateToLocal(pt);
2133     pt->fSort = 0;
2134   }
2135   arr->UnSort();
2136   arr->Sort();  // sorting according z
2137   arr->Expand(arr->GetEntries());
2138   //
2139   //
2140   Int_t nseed=arr->GetEntriesFast();
2141   for (Int_t i=0; i<nseed; i++) {
2142     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2143     if (!pt) continue;
2144     for (Int_t j=0;j<=12;j++){
2145       pt->fOverlapLabels[j] =0;
2146     }
2147   }
2148   for (Int_t i=0; i<nseed; i++) {
2149     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2150     if (!pt) continue;
2151     if (pt->fRemoval>10) continue;
2152     for (Int_t j=i+1; j<nseed; j++){
2153       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2154       //      if (pt2){
2155       if (pt2->fRemoval<=10) {
2156         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2157         SignShared(pt,pt2);
2158       }
2159     }
2160   }
2161 }
2162
2163 void  AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2,  Int_t removalindex)
2164 {
2165   //
2166   //sort trackss according sectors
2167   //
2168   if (fDebug&1) {
2169     Info("RemoveDouble","Number of tracks before double removal- %d\n",arr->GetEntries());
2170   }
2171   //
2172   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2173     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2174     if (!pt) continue;
2175     pt->fSort = 0;
2176   }
2177   arr->UnSort();
2178   arr->Sort();  // sorting according z
2179   arr->Expand(arr->GetEntries());
2180   //
2181   //reset overlap labels
2182   //
2183   Int_t nseed=arr->GetEntriesFast();
2184   for (Int_t i=0; i<nseed; i++) {
2185     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2186     if (!pt) continue;
2187     pt->SetUniqueID(i);
2188     for (Int_t j=0;j<=12;j++){
2189       pt->fOverlapLabels[j] =0;
2190     }
2191   }
2192   //
2193   //sign shared tracks
2194   for (Int_t i=0; i<nseed; i++) {
2195     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2196     if (!pt) continue;
2197     if (pt->fRemoval>10) continue;
2198     Float_t deltac = pt->GetC()*0.1;
2199     for (Int_t j=i+1; j<nseed; j++){
2200       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2201       //      if (pt2){
2202       if (pt2->fRemoval<=10) {
2203         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
2204         if (TMath::Abs(pt->GetC()  -pt2->GetC())>deltac) continue;
2205         if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
2206         //
2207         SignShared(pt,pt2);
2208       }
2209     }
2210   }
2211   //
2212   // remove highly shared tracks
2213   for (Int_t i=0; i<nseed; i++) {
2214     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2215     if (!pt) continue;
2216     if (pt->fRemoval>10) continue;
2217     //
2218     Int_t sumshared =0;
2219     for (Int_t j=0;j<4;j++){
2220       sumshared = pt->fOverlapLabels[j*3+1];
2221     }
2222     Float_t factor = factor1;
2223     if (pt->fRemoval>0) factor = factor2;
2224     if (sumshared/pt->GetNumberOfClusters()>factor){
2225       for (Int_t j=0;j<4;j++){
2226         if (pt->fOverlapLabels[3*j]==0) continue;
2227         if (pt->fOverlapLabels[3*j+1]<5) continue;
2228         if (pt->fRemoval==removalindex) continue;
2229         AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2230         if (!pt2) continue;
2231         if (pt2->GetSigma2C()<pt->GetSigma2C()){
2232           //      pt->fRemoval = removalindex;
2233           delete arr->RemoveAt(i);
2234           break;
2235         }
2236       }
2237     }
2238   }
2239   arr->Compress();
2240   if (fDebug&1) {
2241     Info("RemoveDouble","Number of tracks after double removal- %d\n",arr->GetEntries());
2242   }
2243 }
2244
2245
2246
2247
2248
2249
2250 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2251 {
2252   //
2253   //sort tracks in array according mode criteria
2254   Int_t nseed = arr->GetEntriesFast();
2255   for (Int_t i=0; i<nseed; i++) {
2256     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2257     if (!pt) {
2258       continue;
2259     }
2260     pt->fSort = mode;
2261   }
2262   arr->UnSort();
2263   arr->Sort();
2264 }
2265
2266 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t removalindex)
2267 {
2268
2269   //Loop over all tracks and remove "overlaps"
2270   //
2271   //
2272   Int_t nseed = arr->GetEntriesFast();
2273   Int_t good =0;
2274
2275   for (Int_t i=0; i<nseed; i++) {
2276     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2277     if (!pt) {
2278       delete arr->RemoveAt(i);
2279     }
2280     else{
2281       pt->fSort =1;
2282       pt->fBSigned = kFALSE;
2283     }
2284   }
2285   arr->Compress();
2286   nseed = arr->GetEntriesFast();
2287   arr->UnSort();
2288   arr->Sort();
2289   //
2290   //unsign used
2291   UnsignClusters();
2292   //
2293   for (Int_t i=0; i<nseed; i++) {
2294     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2295     if (!pt) {
2296       continue;
2297     }
2298     Int_t found,foundable,shared;
2299     if (pt->IsActive())
2300       pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2301     else
2302       pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2303     //
2304     Double_t factor = factor2;
2305     if (pt->fBConstrain) factor = factor1;
2306
2307     if ((Float_t(shared)/Float_t(found))>factor){
2308       pt->Desactivate(removalindex);
2309       continue;
2310     }
2311
2312     good++;
2313     for (Int_t i=0; i<160; i++) {
2314       Int_t index=pt->GetClusterIndex2(i);
2315       if (index<0 || index&0x8000 ) continue;
2316       AliTPCclusterMI *c= pt->fClusterPointer[i];
2317       if (!c) continue;
2318       //      if (!c->IsUsed(10)) c->Use(10);
2319       //if (pt->IsActive())
2320       c->Use(10);
2321       //else
2322       //        c->Use(5);
2323     }
2324
2325   }
2326   fNtracks = good;
2327   if (fDebug>0){
2328     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2329   }
2330 }
2331
2332
2333 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2334 {
2335
2336   //Loop over all tracks and remove "overlaps"
2337   //
2338   //
2339   UnsignClusters();
2340   //
2341   Int_t nseed = arr->GetEntriesFast();
2342   Float_t * quality = new Float_t[nseed];
2343   Int_t   * indexes = new Int_t[nseed];
2344   Int_t good =0;
2345   //
2346   //
2347   for (Int_t i=0; i<nseed; i++) {
2348     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2349     if (!pt){
2350       quality[i]=-1;
2351       continue;
2352     }
2353     pt->UpdatePoints();    //select first last max dens points
2354     Float_t * points = pt->GetPoints();
2355     if (points[3]<0.8) quality[i] =-1;
2356     //
2357     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2358   }
2359   TMath::Sort(nseed,quality,indexes);
2360   //
2361   //
2362   for (Int_t itrack=0; itrack<nseed; itrack++) {
2363     Int_t trackindex = indexes[itrack];
2364     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2365     if (quality[trackindex]<0){
2366       if (pt) {
2367         delete arr->RemoveAt(trackindex);
2368       }
2369       else{
2370         arr->RemoveAt(trackindex);
2371       }
2372       continue;
2373     }
2374     //
2375     Int_t first = Int_t(pt->GetPoints()[0]);
2376     Int_t last  = Int_t(pt->GetPoints()[2]);
2377     Double_t factor = (pt->fBConstrain) ? factor1: factor2;
2378     //
2379     Int_t found,foundable,shared;
2380     pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE);
2381     Float_t sharedfactor = Float_t(shared+1)/Float_t(found+1);
2382     //
2383     if (Float_t(shared+1)/Float_t(found+1)>factor){
2384       if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2385       delete arr->RemoveAt(trackindex);
2386       continue;
2387     }
2388
2389     if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
2390       if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2391       delete arr->RemoveAt(trackindex);
2392       continue;
2393     }
2394
2395     good++;
2396     if (sharedfactor>0.4) continue;
2397     for (Int_t i=first; i<last; i++) {
2398       Int_t index=pt->GetClusterIndex2(i);
2399       // if (index<0 || index&0x8000 ) continue;
2400       if (index<0 || index&0x8000 ) continue;
2401       AliTPCclusterMI *c= pt->fClusterPointer[i];
2402       if (!c) continue;
2403       c->Use(10);
2404     }
2405   }
2406   fNtracks = good;
2407   if (fDebug>0){
2408     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2409   }
2410   delete []quality;
2411   delete []indexes;
2412 }
2413
2414 void AliTPCtrackerMI::UnsignClusters()
2415 {
2416   //
2417   // loop over all clusters and unsign them
2418   //
2419
2420   for (Int_t sec=0;sec<fkNIS;sec++){
2421     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2422       AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2423       for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2424         //      if (cl[icl].IsUsed(10))
2425         cl[icl].Use(-1);
2426       cl = fInnerSec[sec][row].fClusters2;
2427       for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2428         //if (cl[icl].IsUsed(10))
2429           cl[icl].Use(-1);
2430     }
2431   }
2432
2433   for (Int_t sec=0;sec<fkNOS;sec++){
2434     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2435       AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2436       for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2437         //if (cl[icl].IsUsed(10))
2438           cl[icl].Use(-1);
2439       cl = fOuterSec[sec][row].fClusters2;
2440       for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2441         //if (cl[icl].IsUsed(10))
2442         cl[icl].Use(-1);
2443     }
2444   }
2445
2446 }
2447
2448
2449
2450 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2451 {
2452   //
2453   //sign clusters to be "used"
2454   //
2455   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2456   // loop over "primaries"
2457
2458   Float_t sumdens=0;
2459   Float_t sumdens2=0;
2460   Float_t sumn   =0;
2461   Float_t sumn2  =0;
2462   Float_t sumchi =0;
2463   Float_t sumchi2 =0;
2464
2465   Float_t sum    =0;
2466
2467   TStopwatch timer;
2468   timer.Start();
2469
2470   Int_t nseed = arr->GetEntriesFast();
2471   for (Int_t i=0; i<nseed; i++) {
2472     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2473     if (!pt) {
2474       continue;
2475     }
2476     if (!(pt->IsActive())) continue;
2477     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2478     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2479       sumdens += dens;
2480       sumdens2+= dens*dens;
2481       sumn    += pt->GetNumberOfClusters();
2482       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2483       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2484       if (chi2>5) chi2=5;
2485       sumchi  +=chi2;
2486       sumchi2 +=chi2*chi2;
2487       sum++;
2488     }
2489   }
2490
2491   Float_t mdensity = 0.9;
2492   Float_t meann    = 130;
2493   Float_t meanchi  = 1;
2494   Float_t sdensity = 0.1;
2495   Float_t smeann    = 10;
2496   Float_t smeanchi  =0.4;
2497
2498
2499   if (sum>20){
2500     mdensity = sumdens/sum;
2501     meann    = sumn/sum;
2502     meanchi  = sumchi/sum;
2503     //
2504     sdensity = sumdens2/sum-mdensity*mdensity;
2505     sdensity = TMath::Sqrt(sdensity);
2506     //
2507     smeann   = sumn2/sum-meann*meann;
2508     smeann   = TMath::Sqrt(smeann);
2509     //
2510     smeanchi = sumchi2/sum - meanchi*meanchi;
2511     smeanchi = TMath::Sqrt(smeanchi);
2512   }
2513
2514
2515   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2516   //
2517   for (Int_t i=0; i<nseed; i++) {
2518     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2519     if (!pt) {
2520       continue;
2521     }
2522     if (pt->fBSigned) continue;
2523     if (pt->fBConstrain) continue;
2524     //if (!(pt->IsActive())) continue;
2525     /*
2526     Int_t found,foundable,shared;
2527     pt->GetClusterStatistic(0,160,found, foundable,shared);
2528     if (shared/float(found)>0.3) {
2529       if (shared/float(found)>0.9 ){
2530         //delete arr->RemoveAt(i);
2531       }
2532       continue;
2533     }
2534     */
2535     Bool_t isok =kFALSE;
2536     if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2537       isok = kTRUE;
2538     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2539       isok =kTRUE;
2540     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2541       isok =kTRUE;
2542     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2543       isok =kTRUE;
2544
2545     if (isok)
2546       for (Int_t i=0; i<160; i++) {
2547         Int_t index=pt->GetClusterIndex2(i);
2548         if (index<0) continue;
2549         AliTPCclusterMI *c= pt->fClusterPointer[i];
2550         if (!c) continue;
2551         //if (!(c->IsUsed(10))) c->Use();
2552         c->Use(10);
2553       }
2554   }
2555
2556
2557   //
2558   Double_t maxchi  = meanchi+2.*smeanchi;
2559
2560   for (Int_t i=0; i<nseed; i++) {
2561     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2562     if (!pt) {
2563       continue;
2564     }
2565     //if (!(pt->IsActive())) continue;
2566     if (pt->fBSigned) continue;
2567     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2568     if (chi>maxchi) continue;
2569
2570     Float_t bfactor=1;
2571     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2572
2573     //sign only tracks with enoug big density at the beginning
2574
2575     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2576
2577
2578     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2579     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2580
2581     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2582     if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2583       minn=0;
2584
2585     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2586       //Int_t noc=pt->GetNumberOfClusters();
2587       pt->fBSigned = kTRUE;
2588       for (Int_t i=0; i<160; i++) {
2589
2590         Int_t index=pt->GetClusterIndex2(i);
2591         if (index<0) continue;
2592         AliTPCclusterMI *c= pt->fClusterPointer[i];
2593         if (!c) continue;
2594         //      if (!(c->IsUsed(10))) c->Use();
2595         c->Use(10);
2596       }
2597     }
2598   }
2599   //  gLastCheck = nseed;
2600   //  arr->Compress();
2601   if (fDebug>0){
2602     timer.Print();
2603   }
2604 }
2605
2606
2607 void  AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2608 {
2609   // stop not active tracks
2610   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2611   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2612   Int_t nseed = arr->GetEntriesFast();
2613   //
2614   for (Int_t i=0; i<nseed; i++) {
2615     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2616     if (!pt) {
2617       continue;
2618     }
2619     if (!(pt->IsActive())) continue;
2620     StopNotActive(pt,row0,th0, th1,th2);
2621   }
2622 }
2623
2624
2625
2626 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2627  Float_t th2) const
2628 {
2629   // stop not active tracks
2630   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2631   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2632   Int_t sumgood1  = 0;
2633   Int_t sumgood2  = 0;
2634   Int_t foundable = 0;
2635   Int_t maxindex = seed->fLastPoint;  //last foundable row
2636   if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2637     seed->Desactivate(10) ;
2638     return;
2639   }
2640
2641   for (Int_t i=row0; i<maxindex; i++){
2642     Int_t index = seed->GetClusterIndex2(i);
2643     if (index!=-1) foundable++;
2644     //if (!c) continue;
2645     if (foundable<=30) sumgood1++;
2646     if (foundable<=50) {
2647       sumgood2++;
2648     }
2649     else{
2650       break;
2651     }
2652   }
2653   if (foundable>=30.){
2654      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2655   }
2656   if (foundable>=50)
2657     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2658 }
2659
2660
2661 Int_t AliTPCtrackerMI::RefitInward(AliESD *event)
2662 {
2663   //
2664   // back propagation of ESD tracks
2665   //
2666   //return 0;
2667   fEvent = event;
2668   ReadSeeds(event,2);
2669   fIteration=2;
2670   //PrepareForProlongation(fSeeds,1);
2671   PropagateForward2(fSeeds);
2672   Int_t ntracks=0;
2673   Int_t nseed = fSeeds->GetEntriesFast();
2674   for (Int_t i=0;i<nseed;i++){
2675     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2676     if (!seed) continue;
2677     if (seed->GetKinkIndex(0)>0)  UpdateKinkQualityD(seed);  // update quality informations for kinks
2678
2679     seed->PropagateTo(fParam->GetInnerRadiusLow());
2680     seed->UpdatePoints();
2681     AliESDtrack *esd=event->GetTrack(i);
2682     seed->CookdEdx(0.02,0.6);
2683     CookLabel(seed,0.1); //For comparison only
2684     if (seed->GetNumberOfClusters()>15){
2685       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2686       esd->SetTPCPoints(seed->GetPoints());
2687       ntracks++;
2688     }
2689     else{
2690       //printf("problem\n");
2691     }
2692   }
2693   //FindKinks(fSeeds,event);
2694   Info("RefitInward","Number of refitted tracks %d",ntracks);
2695   fEvent =0;
2696   //WriteTracks();
2697   return 0;
2698 }
2699
2700
2701 Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2702 {
2703   //
2704   // back propagation of ESD tracks
2705   //
2706
2707   fEvent = event;
2708   fIteration = 1;
2709   ReadSeeds(event,0);
2710   PropagateBack(fSeeds);
2711   Int_t nseed = fSeeds->GetEntriesFast();
2712   Int_t ntracks=0;
2713   for (Int_t i=0;i<nseed;i++){
2714     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2715     if (!seed) continue;
2716     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
2717     seed->UpdatePoints();
2718     AliESDtrack *esd=event->GetTrack(i);
2719     seed->CookdEdx(0.02,0.6);
2720     CookLabel(seed,0.1); //For comparison only
2721     if (seed->GetNumberOfClusters()>15){
2722       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2723       esd->SetTPCPoints(seed->GetPoints());
2724       ntracks++;
2725     }
2726   }
2727   //FindKinks(fSeeds,event);
2728   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2729   fEvent =0;
2730   //WriteTracks();
2731   return 0;
2732 }
2733
2734
2735 void AliTPCtrackerMI::DeleteSeeds()
2736 {
2737   //
2738   //delete Seeds
2739   Int_t nseed = fSeeds->GetEntriesFast();
2740   for (Int_t i=0;i<nseed;i++){
2741     AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2742     if (seed) delete fSeeds->RemoveAt(i);
2743   }
2744   delete fSeeds;
2745   fSeeds =0;
2746 }
2747
2748 void AliTPCtrackerMI::ReadSeeds(AliESD *event, Int_t direction)
2749 {
2750   //
2751   //read seeds from the event
2752
2753   Int_t nentr=event->GetNumberOfTracks();
2754   if (fDebug>0){
2755     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2756   }
2757   if (fSeeds)
2758     DeleteSeeds();
2759   if (!fSeeds){
2760     fSeeds = new TObjArray(nentr);
2761   }
2762   UnsignClusters();
2763   //  Int_t ntrk=0;
2764   for (Int_t i=0; i<nentr; i++) {
2765     AliESDtrack *esd=event->GetTrack(i);
2766     ULong_t status=esd->GetStatus();
2767     if (!(status&AliESDtrack::kTPCin)) continue;
2768     AliTPCtrack t(*esd);
2769     //    AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2770     AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2771     for (Int_t ikink=0;ikink<3;ikink++) {
2772       Int_t index = esd->GetKinkIndex(ikink);
2773       seed->GetKinkIndexes()[ikink] = index;
2774       if (index==0) continue;
2775       index = TMath::Abs(index);
2776       AliESDkink * kink = fEvent->GetKink(index-1);
2777       if (kink&&esd->GetKinkIndex(ikink)<0){
2778         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2779         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,0);
2780       }
2781       if (kink&&esd->GetKinkIndex(ikink)>0){
2782         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2783         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,4);
2784       }
2785
2786     }
2787     if ((status==AliESDtrack::kTPCin)&&(direction==1)) seed->ResetCovariance();
2788     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
2789     if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2790       fSeeds->AddAt(0,i);
2791       delete seed;
2792       continue;
2793     }
2794     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 )  {
2795       Double_t par0[5],par1[5],x;
2796       esd->GetInnerExternalParameters(x,par0);
2797       esd->GetExternalParameters(x,par1);
2798       Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2799       Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2800       Double_t trdchi2=0;
2801       if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2802       //reset covariance if suspicious
2803       if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2804         seed->ResetCovariance();
2805     }
2806
2807     //
2808     //
2809     // rotate to the local coordinate system
2810     //
2811     fSectors=fInnerSec; fN=fkNIS;
2812     Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2813     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2814     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2815     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2816     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2817     if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2818     if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2819     alpha-=seed->GetAlpha();
2820     if (!seed->Rotate(alpha)) {
2821       delete seed;
2822       continue;
2823     }
2824     seed->fEsd = esd;
2825     // sign clusters
2826     for (Int_t irow=0;irow<160;irow++){
2827       Int_t index = seed->GetClusterIndex2(irow);
2828       if (index>0){
2829         //
2830         AliTPCclusterMI * cl = GetClusterMI(index);
2831         seed->fClusterPointer[irow] = cl;
2832         if (cl){
2833           if ((index & 0x8000)==0){
2834             cl->Use(10);  // accepted cluster
2835           }else{
2836             cl->Use(6);   // close cluster not accepted
2837           }
2838         }else{
2839            Info("ReadSeeds","Not found cluster");
2840         }
2841       }
2842     }
2843     fSeeds->AddAt(seed,i);
2844   }
2845 }
2846
2847
2848
2849 //_____________________________________________________________________________
2850 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
2851                                  Float_t deltay, Int_t ddsec) {
2852   //-----------------------------------------------------------------
2853   // This function creates track seeds.
2854   // SEEDING WITH VERTEX CONSTRAIN
2855   //-----------------------------------------------------------------
2856   // cuts[0]   - fP4 cut
2857   // cuts[1]   - tan(phi)  cut
2858   // cuts[2]   - zvertex cut
2859   // cuts[3]   - fP3 cut
2860   Int_t nin0  = 0;
2861   Int_t nin1  = 0;
2862   Int_t nin2  = 0;
2863   Int_t nin   = 0;
2864   Int_t nout1 = 0;
2865   Int_t nout2 = 0;
2866
2867   Double_t x[5], c[15];
2868   //  Int_t di = i1-i2;
2869   //
2870   AliTPCseed * seed = new AliTPCseed;
2871   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2872   Double_t cs=cos(alpha), sn=sin(alpha);
2873   //
2874   //  Double_t x1 =fOuterSec->GetX(i1);
2875   //Double_t xx2=fOuterSec->GetX(i2);
2876
2877   Double_t x1 =GetXrow(i1);
2878   Double_t xx2=GetXrow(i2);
2879
2880   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2881
2882   Int_t imiddle = (i2+i1)/2;    //middle pad row index
2883   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2884   const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2885   //
2886   Int_t ns =sec;
2887
2888   const AliTPCRow& kr1=GetRow(ns,i1);
2889   Double_t ymax  = GetMaxY(i1)-kr1.fDeadZone-1.5;
2890   Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2891
2892   //
2893   // change cut on curvature if it can't reach this layer
2894   // maximal curvature set to reach it
2895   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2896   if (dvertexmax*0.5*cuts[0]>0.85){
2897     cuts[0] = 0.85/(dvertexmax*0.5+1.);
2898   }
2899   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
2900
2901   //  Int_t ddsec = 1;
2902   if (deltay>0) ddsec = 0;
2903   // loop over clusters
2904   for (Int_t is=0; is < kr1; is++) {
2905     //
2906     if (kr1[is]->IsUsed(10)) continue;
2907     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2908     //if (TMath::Abs(y1)>ymax) continue;
2909
2910     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2911
2912     // find possible directions
2913     Float_t anglez = (z1-z3)/(x1-x3);
2914     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
2915     //
2916     //
2917     //find   rotation angles relative to line given by vertex and point 1
2918     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2919     Double_t dvertex  = TMath::Sqrt(dvertex2);
2920     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
2921     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);
2922
2923     //
2924     // loop over 2 sectors
2925     Int_t dsec1=-ddsec;
2926     Int_t dsec2= ddsec;
2927     if (y1<0)  dsec2= 0;
2928     if (y1>0)  dsec1= 0;
2929
2930     Double_t dddz1=0;  // direction of delta inclination in z axis
2931     Double_t dddz2=0;
2932     if ( (z1-z3)>0)
2933       dddz1 =1;
2934     else
2935       dddz2 =1;
2936     //
2937     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2938       Int_t sec2 = sec + dsec;
2939       //
2940       //      AliTPCRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2941       //AliTPCRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2942       AliTPCRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
2943       AliTPCRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2944       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2945       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2946
2947       // rotation angles to p1-p3
2948       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2949       Double_t x2,   y2,   z2;
2950       //
2951       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2952
2953       //
2954       Double_t dxx0 =  (xx2-x3)*cs13r;
2955       Double_t dyy0 =  (xx2-x3)*sn13r;
2956       for (Int_t js=index1; js < index2; js++) {
2957         const AliTPCclusterMI *kcl = kr2[js];
2958         if (kcl->IsUsed(10)) continue;
2959         //
2960         //calcutate parameters
2961         //
2962         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
2963         // stright track
2964         if (TMath::Abs(yy0)<0.000001) continue;
2965         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
2966         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2967         Double_t r02 = (0.25+y0*y0)*dvertex2;
2968         //curvature (radius) cut
2969         if (r02<r2min) continue;
2970
2971         nin0++;
2972         //
2973         Double_t c0  = 1/TMath::Sqrt(r02);
2974         if (yy0>0) c0*=-1.;
2975
2976
2977         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
2978         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2979         Double_t dfi0   = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2980         Double_t dfi1   = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2981         //
2982         //
2983         Double_t z0  =  kcl->GetZ();
2984         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
2985         if (TMath::Abs(zzzz2-z0)>0.5) continue;
2986         nin1++;
2987         //
2988         Double_t dip    = (z1-z0)*c0/dfi1;
2989         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2990         //
2991         y2 = kcl->GetY();
2992         if (dsec==0){
2993           x2 = xx2;
2994           z2 = kcl->GetZ();
2995         }
2996         else
2997           {
2998             // rotation
2999             z2 = kcl->GetZ();
3000             x2= xx2*cs-y2*sn*dsec;
3001             y2=+xx2*sn*dsec+y2*cs;
3002           }
3003
3004         x[0] = y1;
3005         x[1] = z1;
3006         x[2] = x0;
3007         x[3] = dip;
3008         x[4] = c0;
3009         //
3010         //
3011         // do we have cluster at the middle ?
3012         Double_t ym,zm;
3013         GetProlongation(x1,xm,x,ym,zm);
3014         UInt_t dummy;
3015         AliTPCclusterMI * cm=0;
3016         if (TMath::Abs(ym)-ymaxm<0){
3017           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3018           if ((!cm) || (cm->IsUsed(10))) {
3019             continue;
3020           }
3021         }
3022         else{
3023           // rotate y1 to system 0
3024           // get state vector in rotated system
3025           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
3026           Double_t xr2  =  x0*cs+yr1*sn*dsec;
3027           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3028           //
3029           GetProlongation(xx2,xm,xr,ym,zm);
3030           if (TMath::Abs(ym)-ymaxm<0){
3031             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3032             if ((!cm) || (cm->IsUsed(10))) {
3033               continue;
3034             }
3035           }
3036         }
3037
3038
3039         Double_t dym = 0;
3040         Double_t dzm = 0;
3041         if (cm){
3042           dym = ym - cm->GetY();
3043           dzm = zm - cm->GetZ();
3044         }
3045         nin2++;
3046
3047
3048         //
3049         //
3050         Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3051         Double_t sy2=kcl->GetSigmaY2()*2.,     sz2=kcl->GetSigmaZ2()*2.;
3052         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3053         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3054         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3055
3056         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3057         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3058         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3059         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3060         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3061         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3062
3063         Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3064         Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3065         Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3066         Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3067
3068         c[0]=sy1;
3069         c[1]=0.;       c[2]=sz1;
3070         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3071         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3072                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3073         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3074         c[13]=f30*sy1*f40+f32*sy2*f42;
3075         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3076
3077         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3078
3079         UInt_t index=kr1.GetIndex(is);
3080         AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
3081
3082         track->fIsSeeding = kTRUE;
3083         track->fSeed1 = i1;
3084         track->fSeed2 = i2;
3085         track->fSeedType=3;
3086
3087
3088         //if (dsec==0) {
3089           FollowProlongation(*track, (i1+i2)/2,1);
3090           Int_t foundable,found,shared;
3091           track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3092           if ((found<0.55*foundable)  || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3093             seed->Reset();
3094             seed->~AliTPCseed();
3095             continue;
3096           }
3097           //}
3098
3099         nin++;
3100         FollowProlongation(*track, i2,1);
3101
3102
3103         //Int_t rc = 1;
3104         track->fBConstrain =1;
3105         //      track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3106         track->fLastPoint = i1;  // first cluster in track position
3107         track->fFirstPoint = track->fLastPoint;
3108
3109         if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3110             track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3111             track->fNShared>0.4*track->GetNumberOfClusters() ) {
3112           seed->Reset();
3113           seed->~AliTPCseed();
3114           continue;
3115         }
3116         nout1++;
3117         // Z VERTEX CONDITION
3118         Double_t zv;
3119         zv = track->GetZ()+track->GetTgl()/track->GetC()*
3120           ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3121         if (TMath::Abs(zv-z3)>cuts[2]) {
3122           FollowProlongation(*track, TMath::Max(i2-20,0));
3123           zv = track->GetZ()+track->GetTgl()/track->GetC()*
3124             ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3125           if (TMath::Abs(zv-z3)>cuts[2]){
3126             FollowProlongation(*track, TMath::Max(i2-40,0));
3127             zv = track->GetZ()+track->GetTgl()/track->GetC()*
3128               ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
3129             if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
3130               // make seed without constrain
3131               AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3132               FollowProlongation(*track2, i2,1);
3133               track2->fBConstrain = kFALSE;
3134               track2->fSeedType = 1;
3135               arr->AddLast(track2);
3136               seed->Reset();
3137               seed->~AliTPCseed();
3138               continue;
3139             }
3140             else{
3141               seed->Reset();
3142               seed->~AliTPCseed();
3143               continue;
3144
3145             }
3146           }
3147         }
3148
3149         track->fSeedType =0;
3150         arr->AddLast(track);
3151         seed = new AliTPCseed;
3152         nout2++;
3153         // don't consider other combinations
3154         if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
3155           break;
3156       }
3157     }
3158   }
3159   if (fDebug>3){
3160     Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3161   }
3162   delete seed;
3163 }
3164
3165
3166 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3167                                  Float_t deltay) {
3168
3169
3170
3171   //-----------------------------------------------------------------
3172   // This function creates track seeds.
3173   //-----------------------------------------------------------------
3174   // cuts[0]   - fP4 cut
3175   // cuts[1]   - tan(phi)  cut
3176   // cuts[2]   - zvertex cut
3177   // cuts[3]   - fP3 cut
3178
3179
3180   Int_t nin0  = 0;
3181   Int_t nin1  = 0;
3182   Int_t nin2  = 0;
3183   Int_t nin   = 0;
3184   Int_t nout1 = 0;
3185   Int_t nout2 = 0;
3186   Int_t nout3 =0;
3187   Double_t x[5], c[15];
3188   //
3189   // make temporary seed
3190   AliTPCseed * seed = new AliTPCseed;
3191   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3192   //  Double_t cs=cos(alpha), sn=sin(alpha);
3193   //
3194   //
3195
3196   // first 3 padrows
3197   Double_t x1 = GetXrow(i1-1);
3198   const    AliTPCRow& kr1=GetRow(sec,i1-1);
3199   Double_t y1max  = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
3200   //
3201   Double_t x1p = GetXrow(i1);
3202   const    AliTPCRow& kr1p=GetRow(sec,i1);
3203   //
3204   Double_t x1m = GetXrow(i1-2);
3205   const    AliTPCRow& kr1m=GetRow(sec,i1-2);
3206
3207   //
3208   //last 3 padrow for seeding
3209   AliTPCRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
3210   Double_t    x3   =  GetXrow(i1-7);
3211   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3212   //
3213   AliTPCRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
3214   Double_t    x3p   = GetXrow(i1-6);
3215   //
3216   AliTPCRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
3217   Double_t    x3m   = GetXrow(i1-8);
3218
3219   //
3220   //
3221   // middle padrow
3222   Int_t im = i1-4;                           //middle pad row index
3223   Double_t xm         = GetXrow(im);         // radius of middle pad-row
3224   const AliTPCRow& krm=GetRow(sec,im);   //middle pad -row
3225   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3226   //
3227   //
3228   Double_t deltax  = x1-x3;
3229   Double_t dymax   = deltax*cuts[1];
3230   Double_t dzmax   = deltax*cuts[3];
3231   //
3232   // loop over clusters
3233   for (Int_t is=0; is < kr1; is++) {
3234     //
3235     if (kr1[is]->IsUsed(10)) continue;
3236     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3237     //
3238     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
3239     //
3240     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3241     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3242     //
3243     Double_t y3,   z3;
3244     //
3245     //
3246     UInt_t index;
3247     for (Int_t js=index1; js < index2; js++) {
3248       const AliTPCclusterMI *kcl = kr3[js];
3249       if (kcl->IsUsed(10)) continue;
3250       y3 = kcl->GetY();
3251       // apply angular cuts
3252       if (TMath::Abs(y1-y3)>dymax) continue;
3253       x3 = x3;
3254       z3 = kcl->GetZ();
3255       if (TMath::Abs(z1-z3)>dzmax) continue;
3256       //
3257       Double_t angley = (y1-y3)/(x1-x3);
3258       Double_t anglez = (z1-z3)/(x1-x3);
3259       //
3260       Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3261       Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3262       //
3263       Double_t yyym = angley*(xm-x1)+y1;
3264       Double_t zzzm = anglez*(xm-x1)+z1;
3265
3266       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3267       if (!kcm) continue;
3268       if (kcm->IsUsed(10)) continue;
3269
3270       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3271       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3272       //
3273       //
3274       //
3275       Int_t used  =0;
3276       Int_t found =0;
3277       //
3278       // look around first
3279       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3280                                                       anglez*(x1m-x1)+z1,
3281                                                       erry,errz,index);
3282       //
3283       if (kc1m){
3284         found++;
3285         if (kc1m->IsUsed(10)) used++;
3286       }
3287       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3288                                                       anglez*(x1p-x1)+z1,
3289                                                       erry,errz,index);
3290       //
3291       if (kc1p){
3292         found++;
3293         if (kc1p->IsUsed(10)) used++;
3294       }
3295       if (used>1)  continue;
3296       if (found<1) continue;
3297
3298       //
3299       // look around last
3300       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3301                                                       anglez*(x3m-x3)+z3,
3302                                                       erry,errz,index);
3303       //
3304       if (kc3m){
3305         found++;
3306         if (kc3m->IsUsed(10)) used++;
3307       }
3308       else
3309         continue;
3310       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3311                                                       anglez*(x3p-x3)+z3,
3312                                                       erry,errz,index);
3313       //
3314       if (kc3p){
3315         found++;
3316         if (kc3p->IsUsed(10)) used++;
3317       }
3318       else
3319         continue;
3320       if (used>1)  continue;
3321       if (found<3) continue;
3322       //
3323       Double_t x2,y2,z2;
3324       x2 = xm;
3325       y2 = kcm->GetY();
3326       z2 = kcm->GetZ();
3327       //
3328
3329       x[0]=y1;
3330       x[1]=z1;
3331       x[4]=F1(x1,y1,x2,y2,x3,y3);
3332       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3333       nin0++;
3334       //
3335       x[2]=F2(x1,y1,x2,y2,x3,y3);
3336       nin1++;
3337       //
3338       x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3339       //if (TMath::Abs(x[3]) > cuts[3]) continue;
3340       nin2++;
3341       //
3342       //
3343       Double_t sy1=0.1,  sz1=0.1;
3344       Double_t sy2=0.1,  sz2=0.1;
3345       Double_t sy3=0.1,  sy=0.1, sz=0.1;
3346
3347       Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3348       Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3349       Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3350       Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3351       Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3352       Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3353
3354       Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3355       Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3356       Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3357       Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3358
3359       c[0]=sy1;
3360       c[1]=0.;       c[2]=sz1;
3361       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3362       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3363       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3364       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3365       c[13]=f30*sy1*f40+f32*sy2*f42;
3366       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3367
3368       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3369
3370       UInt_t index=kr1.GetIndex(is);
3371       AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3372
3373       track->fIsSeeding = kTRUE;
3374
3375       nin++;
3376       FollowProlongation(*track, i1-7,1);
3377       if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
3378           track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3379         seed->Reset();
3380         seed->~AliTPCseed();
3381         continue;
3382       }
3383       nout1++;
3384       nout2++;
3385       //Int_t rc = 1;
3386       FollowProlongation(*track, i2,1);
3387       track->fBConstrain =0;
3388       track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3389       track->fFirstPoint = track->fLastPoint;
3390
3391       if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3392           track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
3393           track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3394         seed->Reset();
3395         seed->~AliTPCseed();
3396         continue;
3397       }
3398
3399       {
3400         FollowProlongation(*track, TMath::Max(i2-10,0),1);
3401         AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3402         FollowProlongation(*track2, i2,1);
3403         track2->fBConstrain = kFALSE;
3404         track2->fSeedType = 4;
3405         arr->AddLast(track2);
3406         seed->Reset();
3407         seed->~AliTPCseed();
3408       }
3409
3410
3411       //arr->AddLast(track);
3412       //seed = new AliTPCseed;
3413       nout3++;
3414     }
3415   }
3416
3417   if (fDebug>3){
3418     Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3419   }
3420   delete seed;
3421 }
3422
3423
3424 //_____________________________________________________________________________
3425 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3426                                  Float_t deltay, Bool_t /*bconstrain*/) {
3427   //-----------------------------------------------------------------
3428   // This function creates track seeds - without vertex constraint
3429   //-----------------------------------------------------------------
3430   // cuts[0]   - fP4 cut        - not applied
3431   // cuts[1]   - tan(phi)  cut
3432   // cuts[2]   - zvertex cut    - not applied
3433   // cuts[3]   - fP3 cut
3434   Int_t nin0=0;
3435   Int_t nin1=0;
3436   Int_t nin2=0;
3437   Int_t nin3=0;
3438   //  Int_t nin4=0;
3439   //Int_t nin5=0;
3440
3441
3442
3443   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3444   //  Double_t cs=cos(alpha), sn=sin(alpha);
3445   Int_t row0 = (i1+i2)/2;
3446   Int_t drow = (i1-i2)/2;
3447   const AliTPCRow& kr0=fSectors[sec][row0];
3448   AliTPCRow * kr=0;
3449
3450   AliTPCpolyTrack polytrack;
3451   Int_t nclusters=fSectors[sec][row0];
3452   AliTPCseed * seed = new AliTPCseed;
3453
3454   Int_t sumused=0;
3455   Int_t cused=0;
3456   Int_t cnused=0;
3457   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3458     Int_t nfound =0;
3459     Int_t nfoundable =0;
3460     for (Int_t iter =1; iter<2; iter++){   //iterations
3461       const AliTPCRow& krm=fSectors[sec][row0-iter];
3462       const AliTPCRow& krp=fSectors[sec][row0+iter];
3463       const AliTPCclusterMI * cl= kr0[is];
3464
3465       if (cl->IsUsed(10)) {
3466         cused++;
3467       }
3468       else{
3469         cnused++;
3470       }
3471       Double_t x = kr0.GetX();
3472       // Initialization of the polytrack
3473       nfound =0;
3474       nfoundable =0;
3475       polytrack.Reset();
3476       //
3477       Double_t y0= cl->GetY();
3478       Double_t z0= cl->GetZ();
3479       Float_t erry = 0;
3480       Float_t errz = 0;
3481
3482       Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3483       if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;  // seed only at the edge
3484
3485       erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3486       errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3487       polytrack.AddPoint(x,y0,z0,erry, errz);
3488
3489       sumused=0;
3490       if (cl->IsUsed(10)) sumused++;
3491
3492
3493       Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3494       Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3495       //
3496       x = krm.GetX();
3497       AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3498       if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3499         erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3500         errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3501         if (cl1->IsUsed(10))  sumused++;
3502         polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3503       }
3504       //
3505       x = krp.GetX();
3506       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3507       if (cl2) {
3508         erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3509         errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3510         if (cl2->IsUsed(10)) sumused++;
3511         polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3512       }
3513       //
3514       if (sumused>0) continue;
3515       nin0++;
3516       polytrack.UpdateParameters();
3517       // follow polytrack
3518       roadz = 1.2;
3519       roady = 1.2;
3520       //
3521       Double_t yn,zn;
3522       nfoundable = polytrack.GetN();
3523       nfound     = nfoundable;
3524       //
3525       for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3526         Float_t maxdist = 0.8*(1.+3./(ddrow));
3527         for (Int_t delta = -1;delta<=1;delta+=2){
3528           Int_t row = row0+ddrow*delta;
3529           kr = &(fSectors[sec][row]);
3530           Double_t xn = kr->GetX();
3531           Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3532           polytrack.GetFitPoint(xn,yn,zn);
3533           if (TMath::Abs(yn)>ymax) continue;
3534           nfoundable++;
3535           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3536           if (cln) {
3537             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3538             if (dist<maxdist){
3539               /*
3540               erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3541               errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3542               if (cln->IsUsed(10)) {
3543                 //      printf("used\n");
3544                 sumused++;
3545                 erry*=2;
3546                 errz*=2;
3547               }
3548               */
3549               erry=0.1;
3550               errz=0.1;
3551               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3552               nfound++;
3553             }
3554           }
3555         }
3556         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;
3557         polytrack.UpdateParameters();
3558       }
3559     }
3560     if ( (sumused>3) || (sumused>0.5*nfound))  {
3561       //printf("sumused   %d\n",sumused);
3562       continue;
3563     }
3564     nin1++;
3565     Double_t dy,dz;
3566     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3567     AliTPCpolyTrack track2;
3568
3569     polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3570     if (track2.GetN()<0.5*nfoundable) continue;
3571     nin2++;
3572
3573     if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3574       //
3575       // test seed with and without constrain
3576       for (Int_t constrain=0; constrain<=0;constrain++){
3577         // add polytrack candidate
3578
3579         Double_t x[5], c[15];
3580         Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3581         track2.GetBoundaries(x3,x1);
3582         x2 = (x1+x3)/2.;
3583         track2.GetFitPoint(x1,y1,z1);
3584         track2.GetFitPoint(x2,y2,z2);
3585         track2.GetFitPoint(x3,y3,z3);
3586         //
3587         //is track pointing to the vertex ?
3588         Double_t x0,y0,z0;
3589         x0=0;
3590         polytrack.GetFitPoint(x0,y0,z0);
3591
3592         if (constrain) {
3593           x2 = x3;
3594           y2 = y3;
3595           z2 = z3;
3596
3597           x3 = 0;
3598           y3 = 0;
3599           z3 = 0;
3600         }
3601         x[0]=y1;
3602         x[1]=z1;
3603         x[4]=F1(x1,y1,x2,y2,x3,y3);
3604
3605         //      if (TMath::Abs(x[4]) >= cuts[0]) continue;  //
3606         x[2]=F2(x1,y1,x2,y2,x3,y3);
3607
3608         //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3609         //x[3]=F3(x1,y1,x2,y2,z1,z2);
3610         x[3]=F3n(x1,y1,x3,y3,z1,z3,x[4]);
3611         //if (TMath::Abs(x[3]) > cuts[3]) continue;
3612
3613
3614         Double_t sy =0.1, sz =0.1;
3615         Double_t sy1=0.02, sz1=0.02;
3616         Double_t sy2=0.02, sz2=0.02;
3617         Double_t sy3=0.02;
3618
3619         if (constrain){
3620           sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3621         }
3622
3623         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3624         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3625         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3626         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3627         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3628         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3629
3630         Double_t f30=(F3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3631         Double_t f31=(F3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3632         Double_t f32=(F3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3633         Double_t f34=(F3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3634
3635
3636         c[0]=sy1;
3637         c[1]=0.;       c[2]=sz1;
3638         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3639         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3640         c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3641         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3642         c[13]=f30*sy1*f40+f32*sy2*f42;
3643         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3644
3645         //Int_t row1 = fSectors->GetRowNumber(x1);
3646         Int_t row1 = GetRowNumber(x1);
3647
3648         UInt_t index=0;
3649         //kr0.GetIndex(is);
3650         AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3651         track->fIsSeeding = kTRUE;
3652         Int_t rc=FollowProlongation(*track, i2);
3653         if (constrain) track->fBConstrain =1;
3654         else
3655           track->fBConstrain =0;
3656         track->fLastPoint = row1+fInnerSec->GetNRows();  // first cluster in track position
3657         track->fFirstPoint = track->fLastPoint;
3658
3659         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3660             track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3661             track->fNShared>0.4*track->GetNumberOfClusters()) {
3662           //delete track;
3663           seed->Reset();
3664           seed->~AliTPCseed();
3665         }
3666         else {
3667           arr->AddLast(track);
3668           seed = new AliTPCseed;
3669         }
3670         nin3++;
3671       }
3672     }  // if accepted seed
3673   }
3674   if (fDebug>3){
3675     Info("MakeSeeds2","\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3676   }
3677   delete seed;
3678 }
3679
3680
3681 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3682 {
3683   //
3684   //
3685   //reseed using track points
3686   Int_t p0 = int(r0*track->GetNumberOfClusters());     // point 0
3687   Int_t p1 = int(r1*track->GetNumberOfClusters());
3688   Int_t p2 = int(r2*track->GetNumberOfClusters());   // last point
3689   Int_t pp2=0;
3690   Double_t  x0[3],x1[3],x2[3];
3691   x0[0]=-1;
3692   x0[0]=-1;
3693   x0[0]=-1;
3694
3695   // find track position at given ratio of the length
3696   Int_t  sec0, sec1, sec2;
3697   sec0=0;
3698   sec1=0;
3699   sec2=0;
3700   Int_t index=-1;
3701   Int_t clindex;
3702   for (Int_t i=0;i<160;i++){
3703     if (track->fClusterPointer[i]){
3704       index++;
3705       AliTPCTrackerPoint   *trpoint =track->GetTrackPoint(i);
3706       if ( (index<p0) || x0[0]<0 ){
3707         if (trpoint->GetX()>1){
3708           clindex = track->GetClusterIndex2(i);
3709           if (clindex>0){
3710             x0[0] = trpoint->GetX();
3711             x0[1] = trpoint->GetY();
3712             x0[2] = trpoint->GetZ();
3713             sec0  = ((clindex&0xff000000)>>24)%18;
3714           }
3715         }
3716       }
3717
3718       if ( (index<p1) &&(trpoint->GetX()>1)){
3719         clindex = track->GetClusterIndex2(i);
3720         if (clindex>0){
3721           x1[0] = trpoint->GetX();
3722           x1[1] = trpoint->GetY();
3723           x1[2] = trpoint->GetZ();
3724           sec1  = ((clindex&0xff000000)>>24)%18;
3725         }
3726       }
3727       if ( (index<p2) &&(trpoint->GetX()>1)){
3728         clindex = track->GetClusterIndex2(i);
3729         if (clindex>0){
3730           x2[0] = trpoint->GetX();
3731           x2[1] = trpoint->GetY();
3732           x2[2] = trpoint->GetZ();
3733           sec2  = ((clindex&0xff000000)>>24)%18;
3734           pp2 = i;
3735         }
3736       }
3737     }
3738   }
3739
3740   Double_t alpha, cs,sn, xx2,yy2;
3741   //
3742   alpha = (sec1-sec2)*fSectors->GetAlpha();
3743   cs = TMath::Cos(alpha);
3744   sn = TMath::Sin(alpha);
3745   xx2= x1[0]*cs-x1[1]*sn;
3746   yy2= x1[0]*sn+x1[1]*cs;
3747   x1[0] = xx2;
3748   x1[1] = yy2;
3749   //
3750   alpha = (sec0-sec2)*fSectors->GetAlpha();
3751   cs = TMath::Cos(alpha);
3752   sn = TMath::Sin(alpha);
3753   xx2= x0[0]*cs-x0[1]*sn;
3754   yy2= x0[0]*sn+x0[1]*cs;
3755   x0[0] = xx2;
3756   x0[1] = yy2;
3757   //
3758   //
3759   //
3760   Double_t x[5],c[15];
3761   //
3762   x[0]=x2[1];
3763   x[1]=x2[2];
3764   x[4]=F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3765   //  if (x[4]>1) return 0;
3766   x[2]=F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3767   x[3]=F3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3768   //if (TMath::Abs(x[3]) > 2.2)  return 0;
3769   //if (TMath::Abs(x[2]) > 1.99) return 0;
3770   //
3771   Double_t sy =0.1,  sz =0.1;
3772   //
3773   Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3774   Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3775   Double_t sy3=0.01+track->GetSigmaY2();
3776   //
3777   Double_t f40=(F1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3778   Double_t f42=(F1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3779   Double_t f43=(F1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3780   Double_t f20=(F2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3781   Double_t f22=(F2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3782   Double_t f23=(F2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3783   //
3784   Double_t f30=(F3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3785   Double_t f31=(F3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3786   Double_t f32=(F3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3787   Double_t f34=(F3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3788
3789
3790   c[0]=sy1;
3791   c[1]=0.;       c[2]=sz1;
3792   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3793   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3794   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3795   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3796   c[13]=f30*sy1*f40+f32*sy2*f42;
3797   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3798
3799   //  Int_t row1 = fSectors->GetRowNumber(x2[0]);
3800   AliTPCseed *seed=new  AliTPCseed(0, x, c, x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3801   //  Double_t y0,z0,y1,z1, y2,z2;
3802   //seed->GetProlongation(x0[0],y0,z0);
3803   // seed->GetProlongation(x1[0],y1,z1);
3804   //seed->GetProlongation(x2[0],y2,z2);
3805   //  seed =0;
3806   seed->fLastPoint  = pp2;
3807   seed->fFirstPoint = pp2;
3808
3809
3810   return seed;
3811 }
3812
3813
3814 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3815 {
3816   //
3817   //
3818   //reseed using founded clusters
3819   //
3820   // Find the number of clusters
3821   Int_t nclusters = 0;
3822   for (Int_t irow=0;irow<160;irow++){
3823     if (track->GetClusterIndex(irow)>0) nclusters++;
3824   }
3825   //
3826   Int_t ipos[3];
3827   ipos[0] = TMath::Max(int(r0*nclusters),0);             // point 0 cluster
3828   ipos[1] = TMath::Min(int(r1*nclusters),nclusters-1);   //
3829   ipos[2] = TMath::Min(int(r2*nclusters),nclusters-1);   // last point
3830   //
3831   //
3832   Double_t  xyz[3][3];
3833   Int_t     row[3],sec[3]={0,0,0};
3834   //
3835   // find track row position at given ratio of the length
3836   Int_t index=-1;
3837   for (Int_t irow=0;irow<160;irow++){
3838     if (track->GetClusterIndex2(irow)<0) continue;
3839     index++;
3840     for (Int_t ipoint=0;ipoint<3;ipoint++){
3841       if (index<=ipos[ipoint]) row[ipoint] = irow;
3842     }
3843   }
3844   //
3845   //Get cluster and sector position
3846   for (Int_t ipoint=0;ipoint<3;ipoint++){
3847     Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3848     AliTPCclusterMI * cl = GetClusterMI(clindex);
3849     if (cl==0) {
3850       //Error("Bug\n");
3851       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
3852       return 0;
3853     }
3854     sec[ipoint]     = ((clindex&0xff000000)>>24)%18;
3855     xyz[ipoint][0]  = GetXrow(row[ipoint]);
3856     xyz[ipoint][1]  = cl->GetY();
3857     xyz[ipoint][2]  = cl->GetZ();
3858   }
3859   //
3860   //
3861   // Calculate seed state vector and covariance matrix
3862
3863   Double_t alpha, cs,sn, xx2,yy2;
3864   //
3865   alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
3866   cs = TMath::Cos(alpha);
3867   sn = TMath::Sin(alpha);
3868   xx2= xyz[1][0]*cs-xyz[1][1]*sn;
3869   yy2= xyz[1][0]*sn+xyz[1][1]*cs;
3870   xyz[1][0] = xx2;
3871   xyz[1][1] = yy2;
3872   //
3873   alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
3874   cs = TMath::Cos(alpha);
3875   sn = TMath::Sin(alpha);
3876   xx2= xyz[0][0]*cs-xyz[0][1]*sn;
3877   yy2= xyz[0][0]*sn+xyz[0][1]*cs;
3878   xyz[0][0] = xx2;
3879   xyz[0][1] = yy2;
3880   //
3881   //
3882   //
3883   Double_t x[5],c[15];
3884   //
3885   x[0]=xyz[2][1];
3886   x[1]=xyz[2][2];
3887   x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3888   x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
3889   x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
3890   //
3891   Double_t sy =0.1,  sz =0.1;
3892   //
3893   Double_t sy1=0.2, sz1=0.2;
3894   Double_t sy2=0.2, sz2=0.2;
3895   Double_t sy3=0.2;
3896   //
3897   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;
3898   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;
3899   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;
3900   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;
3901   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;
3902   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;
3903   //
3904   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;
3905   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;
3906   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;
3907   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;
3908
3909
3910   c[0]=sy1;
3911   c[1]=0.;       c[2]=sz1;
3912   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3913   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3914   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3915   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3916   c[13]=f30*sy1*f40+f32*sy2*f42;
3917   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3918
3919   //  Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
3920   AliTPCseed *seed=new  AliTPCseed(0, x, c, xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3921   seed->fLastPoint  = row[2];
3922   seed->fFirstPoint = row[2];
3923   return seed;
3924 }
3925
3926
3927 AliTPCseed *AliTPCtrackerMI::ReSeed(AliTPCseed *track,Int_t r0, Bool_t forward)
3928 {
3929   //
3930   //
3931   //reseed using founded clusters
3932   //
3933   Double_t  xyz[3][3];
3934   Int_t     row[3]={0,0,0},sec[3]={0,0,0};
3935   //
3936   // forward direction
3937   if (forward){
3938     for (Int_t irow=r0;irow<160;irow++){
3939       if (track->GetClusterIndex(irow)>0){
3940         row[0] = irow;
3941         break;
3942       }
3943     }
3944     for (Int_t irow=160;irow>r0;irow--){
3945       if (track->GetClusterIndex(irow)>0){
3946         row[2] = irow;
3947         break;
3948       }
3949     }
3950     for (Int_t irow=row[2]-15;irow>row[0];irow--){
3951       if (track->GetClusterIndex(irow)>0){
3952         row[1] = irow;
3953         break;
3954       }
3955     }
3956     //
3957   }
3958   if (!forward){
3959     for (Int_t irow=0;irow<r0;irow++){
3960       if (track->GetClusterIndex(irow)>0){
3961         row[0] = irow;
3962         break;
3963       }
3964     }
3965     for (Int_t irow=r0;irow>0;irow--){
3966       if (track->GetClusterIndex(irow)>0){
3967         row[2] = irow;
3968         break;
3969       }
3970     }
3971     for (Int_t irow=row[2]-15;irow>row[0];irow--){
3972       if (track->GetClusterIndex(irow)>0){
3973         row[1] = irow;
3974         break;
3975       }
3976     }
3977   }
3978   //
3979   if ((row[2]-row[0])<20) return 0;
3980   if (row[1]==0) return 0;
3981   //
3982   //
3983   //Get cluster and sector position
3984   for (Int_t ipoint=0;ipoint<3;ipoint++){
3985     Int_t clindex = track->GetClusterIndex2(row[ipoint]);
3986     AliTPCclusterMI * cl = GetClusterMI(clindex);
3987     if (cl==0) {
3988       //Error("Bug\n");
3989       //      AliTPCclusterMI * cl = GetClusterMI(clindex);
3990       return 0;
3991     }
3992     sec[ipoint]     = ((clindex&0xff000000)>>24)%18;
3993     xyz[ipoint][0]  = GetXrow(row[ipoint]);
3994     AliTPCTrackerPoint * point = track->GetTrackPoint(row[ipoint]);
3995     if (point&&ipoint<2){
3996       //
3997        xyz[ipoint][1]  = point->GetY();
3998        xyz[ipoint][2]  = point->GetZ();
3999     }
4000     else{
4001       xyz[ipoint][1]  = cl->GetY();
4002       xyz[ipoint][2]  = cl->GetZ();
4003     }
4004   }
4005   //
4006   //
4007   //
4008   //
4009   // Calculate seed state vector and covariance matrix
4010
4011   Double_t alpha, cs,sn, xx2,yy2;
4012   //
4013   alpha = (sec[1]-sec[2])*fSectors->GetAlpha();
4014   cs = TMath::Cos(alpha);
4015   sn = TMath::Sin(alpha);
4016   xx2= xyz[1][0]*cs-xyz[1][1]*sn;
4017   yy2= xyz[1][0]*sn+xyz[1][1]*cs;
4018   xyz[1][0] = xx2;
4019   xyz[1][1] = yy2;
4020   //
4021   alpha = (sec[0]-sec[2])*fSectors->GetAlpha();
4022   cs = TMath::Cos(alpha);
4023   sn = TMath::Sin(alpha);
4024   xx2= xyz[0][0]*cs-xyz[0][1]*sn;
4025   yy2= xyz[0][0]*sn+xyz[0][1]*cs;
4026   xyz[0][0] = xx2;
4027   xyz[0][1] = yy2;
4028   //
4029   //
4030   //
4031   Double_t x[5],c[15];
4032   //
4033   x[0]=xyz[2][1];
4034   x[1]=xyz[2][2];
4035   x[4]=F1(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4036   x[2]=F2(xyz[2][0],xyz[2][1],xyz[1][0],xyz[1][1],xyz[0][0],xyz[0][1]);
4037   x[3]=F3n(xyz[2][0],xyz[2][1],xyz[0][0],xyz[0][1],xyz[2][2],xyz[0][2],x[4]);
4038   //
4039   Double_t sy =0.1,  sz =0.1;
4040   //
4041   Double_t sy1=0.2, sz1=0.2;
4042   Double_t sy2=0.2, sz2=0.2;
4043   Double_t sy3=0.2;
4044   //
4045   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;
4046   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;
4047   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;
4048   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;
4049   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;
4050   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;
4051   //
4052   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;
4053   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;
4054   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;
4055   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;
4056
4057
4058   c[0]=sy1;
4059   c[1]=0.;       c[2]=sz1;
4060   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
4061   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
4062   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
4063   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
4064   c[13]=f30*sy1*f40+f32*sy2*f42;
4065   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
4066
4067   //  Int_t row1 = fSectors->GetRowNumber(xyz[2][0]);
4068   AliTPCseed *seed=new  AliTPCseed(0, x, c, xyz[2][0], sec[2]*fSectors->GetAlpha()+fSectors->GetAlphaShift());
4069   seed->fLastPoint  = row[2];
4070   seed->fFirstPoint = row[2];
4071   for (Int_t i=row[0];i<row[2];i++){
4072     seed->fIndex[i] = track->fIndex[i];
4073   }
4074
4075   return seed;
4076 }
4077
4078 void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
4079 {
4080   //
4081   //  find kinks
4082   //
4083   //
4084
4085   TObjArray *kinks= new TObjArray(10000);
4086   Int_t nentries = array->GetEntriesFast();
4087   AliHelix *helixes      = new AliHelix[nentries];
4088   Int_t    *sign         = new Int_t[nentries];
4089   Int_t    *nclusters    = new Int_t[nentries];
4090   Float_t  *alpha        = new Float_t[nentries];
4091   AliESDkink * kink      = new AliESDkink();
4092   Int_t      * usage     = new Int_t[nentries];
4093   Float_t  *zm           = new Float_t[nentries];
4094   Float_t  *z0           = new Float_t[nentries];
4095   Float_t  *fim          = new Float_t[nentries];
4096   Float_t  *shared       = new Float_t[nentries];
4097   Bool_t   *circular     = new Bool_t[nentries];
4098   //
4099   //  nentries = array->GetEntriesFast();
4100   //
4101
4102   //
4103   //
4104   for (Int_t i=0;i<nentries;i++){
4105     sign[i]=0;
4106     usage[i]=0;
4107     AliTPCseed* track = (AliTPCseed*)array->At(i);
4108     if (!track) continue;
4109     shared[i] = kFALSE;
4110     track->UpdatePoints();
4111     if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
4112     }
4113     nclusters[i]=track->GetNumberOfClusters();
4114     alpha[i] = track->GetAlpha();
4115     new (&helixes[i]) AliHelix(*track);
4116     Double_t xyz[3];
4117     helixes[i].Evaluate(0,xyz);
4118     sign[i] = (track->GetC()>0) ? -1:1;
4119     Double_t x,y,z;
4120     x=160;
4121     if (track->GetProlongation(x,y,z)){
4122       zm[i]  = z;
4123       fim[i] = alpha[i]+TMath::ATan2(y,x);
4124     }
4125     else{
4126       zm[i]  = track->GetZ();
4127       fim[i] = alpha[i];
4128     }
4129     z0[i]=1000;
4130     circular[i]= kFALSE;
4131     if (track->GetProlongation(0,y,z))  z0[i] = TMath::Abs(z);
4132   }
4133   //
4134   //
4135   TStopwatch timer;
4136   timer.Start();
4137   Int_t ncandidates =0;
4138   Int_t nall =0;
4139   Int_t ntracks=0;
4140   Double_t phase[2][2],radius[2];
4141
4142   //
4143   // Find circling track
4144   //  TTreeSRedirector cstream("circling.root");
4145   //
4146   for (Int_t i0=0;i0<nentries;i0++){
4147     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
4148     if (!track0) continue;
4149     if (track0->fN<40) continue;
4150     if (TMath::Abs(1./track0->fP4)>200) continue;
4151     for (Int_t i1=i0+1;i1<nentries;i1++){
4152       AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
4153       if (!track1) continue;
4154       if (TMath::Abs(1./track1->fP4)>200) continue;
4155       if (track1->fN<40) continue;
4156       if (z0[i0]<20&&z0[i1]<20) continue;
4157       if (track1->fP4*track0->fP4>0) continue;
4158       if (track1->fP3*track0->fP3>0) continue;
4159       if (track0->fBConstrain&&track1->fBConstrain) continue;
4160       if (max(TMath::Abs(1./track0->fP4),TMath::Abs(1./track1->fP4))>190) continue;
4161       if (track0->fBConstrain&&TMath::Abs(track1->fP4)<TMath::Abs(track0->fP4)) continue; //returning - lower momenta
4162       if (track1->fBConstrain&&TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)) continue; //returning - lower momenta
4163       //
4164       if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue;
4165       Float_t xc0 = helixes[i0].GetHelix(6);
4166       Float_t yc0 = helixes[i0].GetHelix(7);
4167       Float_t r0  = helixes[i0].GetHelix(8);
4168       Float_t xc1 = helixes[i1].GetHelix(6);
4169       Float_t yc1 = helixes[i1].GetHelix(7);
4170       Float_t r1  = helixes[i1].GetHelix(8);
4171
4172       Float_t rmean = (r0+r1)*0.5;
4173       Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
4174       if (delta>30) continue;
4175       if (delta>rmean*0.25) continue;
4176       if (TMath::Abs(r0-r1)/rmean>0.3) continue;
4177       //
4178       Int_t npoints = helixes[i0].GetRPHIintersections(helixes[i1], phase, radius,10);
4179       if (npoints==0) continue;
4180       helixes[i0].GetClosestPhases(helixes[i1], phase);
4181       //
4182       Double_t xyz0[3];
4183       Double_t xyz1[3];
4184       Double_t hangles[3];
4185       helixes[i0].Evaluate(phase[0][0],xyz0);
4186       helixes[i1].Evaluate(phase[0][1],xyz1);
4187
4188       helixes[i0].GetAngle(phase[0][0],helixes[i1],phase[0][1],hangles);
4189       Double_t deltah[2],deltabest;
4190       if (hangles[2]<2.8) continue;
4191       /*
4192       cstream<<"C"<<track0->fLab<<track1->fLab<<
4193         track0->fP3<<track1->fP3<<
4194         track0->fP4<<track1->fP4<<
4195         delta<<rmean<<npoints<<
4196         hangles[0]<<hangles[2]<<
4197         xyz0[2]<<xyz1[2]<<radius[0]<<"\n";
4198       */
4199       if (npoints>0){
4200         Int_t ibest=0;
4201         helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],5);
4202         if (npoints==2){
4203           helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],5);
4204           if (deltah[1]<deltah[0]) ibest=1;
4205         }
4206         deltabest  = TMath::Sqrt(deltah[ibest]);
4207         helixes[i0].Evaluate(phase[ibest][0],xyz0);
4208         helixes[i1].Evaluate(phase[ibest][1],xyz1);
4209         helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
4210         //Double_t radiusbest = TMath::Sqrt(radius[ibest]);
4211         //
4212         if (deltabest<25){
4213           /*
4214           cstream<<"CC1"<<track0->fLab<<track1->fLab<<   //0-1
4215             track0->fP3<<track1->fP3<<                   //2-3
4216             track0->fP4<<track1->fP4<<                   //4-5
4217             delta<<rmean<<npoints<<                      //6-8
4218             hangles[0]<<hangles[2]<<                     //9-10
4219             xyz0[2]<<xyz1[2]<<radiusbest<<deltabest<< //11--14
4220             track0->fFirstPoint<<track1->fFirstPoint<<   //15-16
4221             track0->fLastPoint<<track1->fLastPoint<<     //17-18
4222             track0<<track1<<                             //19-20
4223             phase[ibest][0]<<phase[ibest][1]<<"\n"; //21-22
4224           */
4225         }
4226         Bool_t sign =kFALSE;
4227         if (deltabest<3&&hangles[2]>3.06) sign =kTRUE;
4228         if (TMath::Min(TMath::Abs(track0->GetD()),TMath::Abs(track1->GetD()))>10&&deltabest<6&&hangles[2]>3.) sign= kTRUE;
4229         if (sign){
4230           circular[i0] = kTRUE;
4231           circular[i1] = kTRUE;
4232         }
4233       }
4234     }
4235   }
4236   //
4237   //
4238   //
4239   //
4240   for (Int_t i =0;i<nentries;i++){
4241     if (sign[i]==0) continue;
4242     AliTPCseed * track0 = (AliTPCseed*)array->At(i);
4243     ntracks++;
4244     //
4245     Double_t cradius0 = 40*40;
4246     Double_t cradius1 = 270*270;
4247     Double_t cdist1=8.;
4248     Double_t cdist2=8.;
4249     Double_t cdist3=0.55;
4250     for (Int_t j =i+1;j<nentries;j++){
4251       nall++;
4252       if (sign[j]*sign[i]<1) continue;
4253       if ( (nclusters[i]+nclusters[j])>200) continue;
4254       if ( (nclusters[i]+nclusters[j])<80) continue;
4255       if ( TMath::Abs(zm[i]-zm[j])>60.) continue;
4256       if ( TMath::Abs(fim[i]-fim[j])>0.6 && TMath::Abs(fim[i]-fim[j])<5.7 ) continue;
4257       //AliTPCseed * track1 = (AliTPCseed*)array->At(j);  Double_t phase[2][2],radius[2];
4258       Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4259       if (npoints<1) continue;
4260       // cuts on radius
4261       if (npoints==1){
4262         if (radius[0]<cradius0||radius[0]>cradius1) continue;
4263       }
4264       else{
4265         if ( (radius[0]<cradius0||radius[0]>cradius1) && (radius[1]<cradius0||radius[1]>cradius1) ) continue;
4266       }
4267       //
4268       Double_t delta1=10000,delta2=10000;
4269       // cuts on the intersection radius
4270       helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4271       if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4272       if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4273       if (npoints==2){
4274         helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);
4275         if (radius[1]<20&&delta2<1) continue;  //intersection at vertex
4276         if (radius[1]<10&&delta2<3) continue;  //intersection at vertex
4277       }
4278       //
4279       Double_t distance1 = TMath::Min(delta1,delta2);
4280       if (distance1>cdist1) continue;  // cut on DCA linear approximation
4281       //
4282       npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,20);
4283       helixes[i].ParabolicDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta1);
4284       if (radius[0]<20&&delta1<1) continue; //intersection at vertex
4285       if (radius[0]<10&&delta1<3) continue; //intersection at vertex
4286       //
4287       if (npoints==2){
4288         helixes[i].ParabolicDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta2);