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