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 */
  19
  20
  21
  22
  23
  24
  25
  26 /*
  27   AliTPC parallel tracker -
  28   How to use?  -
  29   run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
  30   run AliTPCFindTracksMI.C macro - to find tracks
  31   tracks are written to AliTPCtracks.root file
  32   for comparison also seeds are written to the same file - to special branch
  33 */
  34
  35 //-------------------------------------------------------
  36 //          Implementation of the TPC tracker
  37 //
  38 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
  39 //
  40 //-------------------------------------------------------
  41
  42 #include <TObjArray.h>
  43 #include <TFile.h>
  44 #include <TTree.h>
  45 #include "Riostream.h"
  46
  47 #include "AliTPCtrackerMI.h"
  48 #include "AliTPCclusterMI.h"
  49 #include "AliTPCParam.h"
  50 #include "AliTPCClustersRow.h"
  51 #include "AliComplexCluster.h"
  52 #include "AliTPCpolyTrack.h"
  53 #include "TStopwatch.h"
  54
  55
  56 ClassImp(AliTPCseed)
  57 ClassImp(AliTPCKalmanSegment)
  58
  59
  60
  61 //_____________________________________________________________________________
  62
  63 AliTPCKalmanSegment::AliTPCKalmanSegment(){
  64   //
  65   //
  66   fX=fAlpha=fChi2=0;
  67   for (Int_t i=0;i<5;i++) fState[i] = 0.;
  68   for (Int_t i=0;i<15;i++) fCovariance[i] = 0.;
  69   fNCFoundable = 0;
  70   fNC          = 0;
  71   //  fN           = 0;
  72 }
  73
  74 //_____________________________________________________________________________
  75
  76 void AliTPCKalmanSegment::Init(AliTPCseed* seed)
  77 {
  78   // in initialization
  79   // initial entrance integral chi2, fNCFoundable and fNC stored
  80   fNCFoundable = seed->fNFoundable;
  81   fNC          = seed->GetNumberOfClusters();
  82   fChi2        = seed->GetChi2();
  83 }
  84
  85
  86 void AliTPCKalmanSegment::Finish(AliTPCseed* seed)
  87 {
  88   //
  89   // in finish state vector stored and chi2 and fNC... calculated
  90   Double_t x;
  91   Double_t state[5];
  92   Double_t cov[15];
  93   seed->GetExternalParameters(x,state);
  94   seed->GetExternalCovariance(cov);
  95   //float precision for tracklet
  96   for (Int_t i=0;i<5;i++) fState[i] = state[i];
  97   for (Int_t i=0;i<15;i++) fCovariance[i] = cov[i];
  98   //
  99   // in current   seed integral track characteristic
 100   // for tracklet differenciation between beginning (Init state) and Finish state
 101   fNCFoundable =  seed->fNFoundable - fNCFoundable;
 102   fNC          =  seed->GetNumberOfClusters() - fNC;
 103   fChi2        =  seed->GetChi2()-fChi2;
 104   //
 105 }
 106
 107 void AliTPCKalmanSegment::GetState(Double_t &x, Double_t & alpha, Double_t state[5])
 108 {
 109   //
 110   x = fX;
 111   alpha = fAlpha;
 112   for (Int_t i=0;i<5;i++) state[i] = fState[i];
 113 }
 114
 115 void AliTPCKalmanSegment::GetCovariance(Double_t covariance[15])
 116 {
 117   //
 118   for (Int_t i=0;i<5;i++) covariance[i] = fCovariance[i];
 119 }
 120
 121 void AliTPCKalmanSegment::GetStatistic(Int_t & nclusters, Int_t & nfoundable, Float_t & chi2)
 122 {
 123   //
 124   //
 125   nclusters  = fNC;
 126   nfoundable = fNCFoundable;
 127   chi2       = fChi2;
 128 }
 129
 130
 131
 132 AliTPCclusterTracks::AliTPCclusterTracks(){
 133   // class for storing overlaping info
 134   fTrackIndex[0]=-1;
 135   fTrackIndex[1]=-1;
 136   fTrackIndex[2]=-1;
 137   fDistance[0]=1000;
 138   fDistance[1]=1000;
 139   fDistance[2]=1000;
 140 }
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){
 151
 152   Int_t sec=(i&0xff000000)>>24;
 153   Int_t row = (i&0x00ff0000)>>16;
 154   track->fRow=(i&0x00ff0000)>>16;
 155   track->fSector = sec;
 156   //  Int_t index = i&0xFFFF;
 157   if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
 158   track->fClusterIndex[track->fRow] = i;
 159   track->fFirstPoint = row;
 160   if ( track->fLastPoint<row) track->fLastPoint =row;
 161   //
 162
 163   AliTPCTrackPoint   *trpoint =track->GetTrackPoint(track->fRow);
 164   Float_t angle2 = track->GetSnp()*track->GetSnp();
 165   angle2 = TMath::Sqrt(angle2/(1-angle2));
 166   //
 167   //SET NEW Track Point
 168   //
 169   if (c!=0){
 170     //if we have a cluster
 171     trpoint->GetCPoint().SetY(c->GetY());
 172     trpoint->GetCPoint().SetZ(c->GetZ());
 173     //
 174     trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY));
 175     trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ));
 176     //
 177     trpoint->GetCPoint().SetType(c->GetType());
 178     trpoint->GetCPoint().SetQ(c->GetQ());
 179     trpoint->GetCPoint().SetMax(c->GetMax());
 180     //
 181     trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2));
 182     trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2));
 183     //
 184   }
 185   trpoint->GetTPoint().SetX(track->GetX());
 186   trpoint->GetTPoint().SetY(track->GetY());
 187   trpoint->GetTPoint().SetZ(track->GetZ());
 188   //
 189   trpoint->GetTPoint().SetAngleY(angle2);
 190   trpoint->GetTPoint().SetAngleZ(track->GetTgl());
 191
 192
 193
 194   if (chi2>10){
 195     //    printf("suspicious chi2 %f\n",chi2);
 196   }
 197   //  if (track->fIsSeeding){
 198     track->fErrorY2 *= 1.2;
 199     track->fErrorY2 += 0.0064;
 200     track->fErrorZ2 *= 1.2;
 201     track->fErrorY2 += 0.005;
 202
 203     //}
 204
 205   return track->Update(c,chi2,i);
 206
 207 }
 208 //_____________________________________________________________________________
 209 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn):
 210 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 211 {
 212   //---------------------------------------------------------------------
 213   // The main TPC tracker constructor
 214   //---------------------------------------------------------------------
 215   fInnerSec=new AliTPCSector[fkNIS];
 216   fOuterSec=new AliTPCSector[fkNOS];
 217
 218   Int_t i;
 219   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
 220   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
 221
 222   fN=0;  fSectors=0;
 223
 224   fClustersArray.Setup(par);
 225   fClustersArray.SetClusterType("AliTPCclusterMI");
 226
 227   char   cname[100];
 228   if (eventn==-1) {
 229     sprintf(cname,"TreeC_TPC");
 230   }
 231   else {
 232     sprintf(cname,"TreeC_TPC_%d",eventn);
 233   }
 234
 235   fClustersArray.ConnectTree(cname);
 236
 237   fEventN = eventn;
 238   fSeeds=0;
 239   fNtracks = 0;
 240   fParam = par;
 241 }
 242
 243 //_____________________________________________________________________________
 244 AliTPCtrackerMI::~AliTPCtrackerMI() {
 245   //------------------------------------------------------------------
 246   // TPC tracker destructor
 247   //------------------------------------------------------------------
 248   delete[] fInnerSec;
 249   delete[] fOuterSec;
 250   if (fSeeds) {
 251     fSeeds->Delete();
 252     delete fSeeds;
 253   }
 254 }
 255
 256
 257 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
 258   //
 259   //
 260   Float_t snoise2;
 261   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 262
 263   //cluster "quality"
 264   Float_t rsigmay = 1;
 265   Int_t ctype = 0;
 266
 267   //standard if we don't have cluster - take MIP
 268   const Float_t chmip = 50.;
 269   Float_t amp = chmip/0.3;
 270   Float_t nel;
 271   Float_t nprim;
 272   if (cl){
 273     amp = cl->GetQ();
 274     rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY);
 275     ctype = cl->GetType();
 276   }
 277
 278
 279   Float_t landau=2 ;    //landau fluctuation part
 280   Float_t gg=2;         // gg fluctuation part
 281   Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
 282
 283
 284   if (fSectors==fInnerSec){
 285     snoise2 = 0.0004;
 286     nel     = 0.268*amp;
 287     nprim   = 0.155*amp;
 288     gg      = (2+0.0002*amp)/nel;
 289     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
 290     if (landau>1) landau=1;
 291   }
 292   else {
 293     snoise2 = 0.0004;
 294     nel     = 0.3*amp;
 295     nprim   = 0.133*amp;
 296     gg      = (2+0.0002*amp)/nel;
 297     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
 298     if (landau>1) landau=1;
 299   }
 300
 301
 302   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
 303   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 304   angle2 = angle2/(1-angle2);
 305   Float_t angular = landau*angle2*padlength*padlength/12.;
 306   Float_t res = sdiff + angular;
 307
 308
 309   if ((ctype==0) && (fSectors ==fOuterSec))
 310     res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2));
 311
 312   if ((ctype==0) && (fSectors ==fInnerSec))
 313     res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2));
 314
 315
 316   if ((ctype>0))
 317     res*= TMath::Power((rsigmay+0.5),1.5)+0.0064;
 318
 319   if (ctype<0)
 320     res*=2.4;  // overestimate error 2 times
 321
 322   res+= snoise2;
 323
 324   if (res<2*snoise2)
 325     res = 2*snoise2;
 326
 327   seed->SetErrorY2(res);
 328   return res;
 329 }
 330
 331
 332
 333
 334 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 335   //
 336   //
 337   Float_t snoise2;
 338   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 339   //signal quality
 340   Float_t rsigmaz=1;
 341   Int_t ctype =0;
 342
 343   const Float_t chmip = 50.;
 344   Float_t amp = chmip/0.3;
 345   Float_t nel;
 346   Float_t nprim;
 347   if (cl){
 348     amp = cl->GetQ();
 349     rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ);
 350     ctype = cl->GetType();
 351   }
 352
 353   //
 354   Float_t landau=2 ;    //landau fluctuation part
 355   Float_t gg=2;         // gg fluctuation part
 356   Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
 357
 358   if (fSectors==fInnerSec){
 359     snoise2 = 0.0004;
 360     nel     = 0.268*amp;
 361     nprim   = 0.155*amp;
 362     gg      = (2+0.0002*amp)/nel;
 363     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
 364     if (landau>1) landau=1;
 365   }
 366   else {
 367     snoise2 = 0.0004;
 368     nel     = 0.3*amp;
 369     nprim   = 0.133*amp;
 370     gg      = (2+0.0002*amp)/nel;
 371     landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
 372     if (landau>1) landau=1;
 373   }
 374   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
 375
 376   Float_t angle = seed->GetTgl();
 377   Float_t angular = landau*angle*angle*padlength*padlength/12.;
 378   Float_t res = sdiff + angular;
 379
 380   if ((ctype==0) && (fSectors ==fOuterSec))
 381     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
 382
 383   if ((ctype==0) && (fSectors ==fInnerSec))
 384     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
 385   if ((ctype>0))
 386     res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064;  //0.31+0.147*ctype;
 387   if (ctype<0)
 388     res*=1.3;
 389   if ((ctype<0) &&amp<70)
 390     res*=1.3;
 391
 392   res += snoise2;
 393   if (res<2*snoise2)
 394      res = 2*snoise2;
 395
 396   seed->SetErrorZ2(res);
 397   return res;
 398 }
 399
 400
 401
 402
 403
 404 void AliTPCseed::Reset()
 405 {
 406   //
 407   //PH  SetN(0);
 408   fNFoundable = 0;
 409   ResetCovariance();
 410   SetChi2(0);
 411   for (Int_t i=0;i<200;i++) fClusterIndex[i]=-1;
 412 }
 413
 414
 415 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
 416 {
 417   //-----------------------------------------------------------------
 418   // This function find proloncation of a track to a reference plane x=xk.
 419   // doesn't change internal state of the track
 420   //-----------------------------------------------------------------
 421
 422   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
 423   //  Double_t y1=fP0, z1=fP1;
 424   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
 425   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
 426
 427   y = fP0;
 428   z = fP1;
 429   y += dx*(c1+c2)/(r1+r2);
 430   z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 431   return 0;
 432 }
 433
 434
 435 //_____________________________________________________________________________
 436 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
 437 {
 438   //-----------------------------------------------------------------
 439   // This function calculates a predicted chi2 increment.
 440   //-----------------------------------------------------------------
 441   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
 442   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
 443   r00+=fC00; r01+=fC10; r11+=fC11;
 444
 445   Double_t det=r00*r11 - r01*r01;
 446   if (TMath::Abs(det) < 1.e-10) {
 447     Int_t n=GetNumberOfClusters();
 448     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
 449     return 1e10;
 450   }
 451   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 452
 453   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
 454
 455   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
 456 }
 457
 458
 459 //_________________________________________________________________________________________
 460
 461
 462 Int_t AliTPCseed::Compare(const TObject *o) const {
 463   //-----------------------------------------------------------------
 464   // This function compares tracks according to the sector - for given sector according z
 465   //-----------------------------------------------------------------
 466   AliTPCseed *t=(AliTPCseed*)o;
 467   if (t->fRelativeSector>fRelativeSector) return -1;
 468   if (t->fRelativeSector<fRelativeSector) return 1;
 469
 470   Double_t z2 = t->GetZ();
 471   Double_t z1 = GetZ();
 472   if (z2>z1) return 1;
 473   if (z2<z1) return -1;
 474   return 0;
 475 }
 476
 477 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
 478 {
 479   //rotate to track "local coordinata
 480   Float_t x = seed->GetX();
 481   Float_t y = seed->GetY();
 482   Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
 483   if (y > ymax) {
 484     seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
 485     if (!seed->Rotate(fSectors->GetAlpha()))
 486       return;
 487   } else if (y <-ymax) {
 488     seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
 489     if (!seed->Rotate(-fSectors->GetAlpha()))
 490       return;
 491   }
 492
 493 }
 494
 495
 496
 497
 498 //_____________________________________________________________________________
 499 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
 500   //-----------------------------------------------------------------
 501   // This function associates a cluster with this track.
 502   //-----------------------------------------------------------------
 503   //  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
 504   //Double_t r00=sigmay2, r01=0., r11=sigmaz2;
 505   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
 506
 507   r00+=fC00; r01+=fC10; r11+=fC11;
 508   Double_t det=r00*r11 - r01*r01;
 509   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 510
 511   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
 512   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
 513   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
 514   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
 515   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
 516
 517   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
 518   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
 519   if (TMath::Abs(cur*fX-eta) >= 0.9) {
 520     //    Int_t n=GetNumberOfClusters();
 521     //if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
 522     return 0;
 523   }
 524
 525   fP0 += k00*dy + k01*dz;
 526   fP1 += k10*dy + k11*dz;
 527   fP2  = eta;
 528   fP3 += k30*dy + k31*dz;
 529   fP4  = cur;
 530
 531   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
 532   Double_t c12=fC21, c13=fC31, c14=fC41;
 533
 534   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
 535   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
 536   fC40-=k00*c04+k01*c14;
 537
 538   fC11-=k10*c01+k11*fC11;
 539   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
 540   fC41-=k10*c04+k11*c14;
 541
 542   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
 543   fC42-=k20*c04+k21*c14;
 544
 545   fC33-=k30*c03+k31*c13;
 546   fC43-=k40*c03+k41*c13;
 547
 548   fC44-=k40*c04+k41*c14;
 549
 550   Int_t n=GetNumberOfClusters();
 551   fIndex[n]=index;
 552   SetNumberOfClusters(n+1);
 553   SetChi2(GetChi2()+chisq);
 554
 555   return 1;
 556 }
 557
 558
 559
 560 //_____________________________________________________________________________
 561 Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
 562                    Double_t x2,Double_t y2,
 563                    Double_t x3,Double_t y3)
 564 {
 565   //-----------------------------------------------------------------
 566   // Initial approximation of the track curvature
 567   //-----------------------------------------------------------------
 568   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 569   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 570                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 571   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 572                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 573
 574   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 575
 576   return -xr*yr/sqrt(xr*xr+yr*yr);
 577 }
 578
 579
 580 //_____________________________________________________________________________
 581 Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
 582                    Double_t x2,Double_t y2,
 583                    Double_t x3,Double_t y3)
 584 {
 585   //-----------------------------------------------------------------
 586   // Initial approximation of the track curvature times center of curvature
 587   //-----------------------------------------------------------------
 588   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 589   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 590                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 591   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 592                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 593
 594   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 595
 596   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 597 }
 598
 599 //_____________________________________________________________________________
 600 Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1,
 601                    Double_t x2,Double_t y2,
 602                    Double_t z1,Double_t z2)
 603 {
 604   //-----------------------------------------------------------------
 605   // Initial approximation of the tangent of the track dip angle
 606   //-----------------------------------------------------------------
 607   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 608 }
 609
 610
 611 void AliTPCtrackerMI::LoadClusters()
 612 {
 613   //
 614   // load clusters to the memory
 615   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 616   for (Int_t i=0; i<j; i++) {
 617     fClustersArray.LoadEntry(i);
 618   }
 619 }
 620
 621 void AliTPCtrackerMI::UnloadClusters()
 622 {
 623   //
 624   // load clusters to the memory
 625   Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 626   for (Int_t i=0; i<j; i++) {
 627     fClustersArray.ClearSegment(i);
 628   }
 629 }
 630
 631
 632
 633 //_____________________________________________________________________________
 634 void AliTPCtrackerMI::LoadOuterSectors() {
 635   //-----------------------------------------------------------------
 636   // This function fills outer TPC sectors with clusters.
 637   //-----------------------------------------------------------------
 638   UInt_t index;
 639   //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 640   Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
 641   for (Int_t i=0; i<j; i++) {
 642     //  AliSegmentID *s=fClustersArray.LoadEntry(i);
 643     AliSegmentID *s= const_cast<AliSegmentID*>(fClustersArray.At(i));
 644     if (!s) continue;
 645     Int_t sec,row;
 646     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 647     par->AdjustSectorRow(s->GetID(),sec,row);
 648     if (sec<fkNIS*2) continue;
 649     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 650     Int_t ncl=clrow->GetArray()->GetEntriesFast();
 651     while (ncl--) {
 652       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
 653       index=(((sec<<8)+row)<<16)+ncl;
 654       fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
 655     }
 656   }
 657   fN=fkNOS;
 658   fSectors=fOuterSec;
 659 }
 660
 661
 662 //_____________________________________________________________________________
 663 void AliTPCtrackerMI::LoadInnerSectors() {
 664   //-----------------------------------------------------------------
 665   // This function fills inner TPC sectors with clusters.
 666   //-----------------------------------------------------------------
 667   UInt_t index;
 668   //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 669   Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
 670   for (Int_t i=0; i<j; i++) {
 671     //   AliSegmentID *s=fClustersArray.LoadEntry(i);
 672     AliSegmentID *s=const_cast<AliSegmentID*>(fClustersArray.At(i));
 673     if (!s) continue;
 674     Int_t sec,row;
 675     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 676     par->AdjustSectorRow(s->GetID(),sec,row);
 677     if (sec>=fkNIS*2) continue;
 678     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 679     Int_t ncl=clrow->GetArray()->GetEntriesFast();
 680     while (ncl--) {
 681       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
 682       index=(((sec<<8)+row)<<16)+ncl;
 683       fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
 684     }
 685   }
 686
 687   fN=fkNIS;
 688   fSectors=fInnerSec;
 689 }
 690
 691 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
 692   //-----------------------------------------------------------------
 693   // This function tries to find a track prolongation to next pad row
 694   //-----------------------------------------------------------------
 695   //  Double_t xt=t.GetX();
 696   //  Int_t row = fSectors->GetRowNumber(xt)-1;
 697   //  if (row < nr) return 1; // don't prolongate if not information until now -
 698   //
 699   Double_t  x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
 700   //  if (t.GetRadius()>x+10 ) return 0;
 701
 702   if (!t.PropagateTo(x)) {
 703     t.fStopped = kTRUE;
 704     return 0;
 705   }
 706   // update current
 707   t.fCurrentSigmaY = GetSigmaY(&t);
 708   t.fCurrentSigmaZ = GetSigmaZ(&t);
 709   //
 710   AliTPCclusterMI *cl=0;
 711   UInt_t index=0;
 712   const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
 713   Double_t sy2=ErrY2(&t)*2;
 714   Double_t sz2=ErrZ2(&t)*2;
 715
 716
 717   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
 718   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
 719   Double_t  y=t.GetY(), z=t.GetZ();
 720
 721   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
 722     t.fInDead = kTRUE;
 723     Int_t row = nr;
 724     if (fSectors==fOuterSec) row += fParam->GetNRowLow();
 725     t.fClusterIndex[row] = -1;
 726     return 0;
 727   }
 728   else
 729     {
 730       if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
 731       else
 732         return 0;
 733     }
 734   //calculate
 735   Float_t maxdistance = roady*roady + roadz*roadz;
 736   if (krow) {
 737     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
 738       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
 739       if (c->GetZ() > z+roadz) break;
 740       if ( (c->GetY()-y) >  roady ) continue;
 741       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
 742       if (maxdistance>distance) {
 743         maxdistance = distance;
 744         cl=c;
 745         //      index=krow.GetIndex(i);
 746         index =i;
 747       }
 748     }
 749   }
 750   if (cl) {
 751     //    Double_t sy2= ErrY2(&t,cl);
 752     //    Double_t sz2= ErrZ2(&t,cl);
 753     //    Double_t chi2= t.GetPredictedChi2(cl);
 754     //    UpdateTrack(&t,cl,chi2,index);
 755
 756     t.fCurrentCluster = cl;
 757     t.fCurrentClusterIndex1 = krow.GetIndex(index);
 758     t.fCurrentClusterIndex2 = index;
 759     Double_t sy2=ErrY2(&t,t.fCurrentCluster);
 760     Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
 761
 762     Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
 763     Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
 764
 765     Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
 766     Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
 767
 768     Double_t rdistance  = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
 769
 770
 771     //    printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
 772     if ( (rdistancey>1) || (rdistancez>1)) return 0;
 773     if (rdistance>4) return 0;
 774
 775     if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)
 776         return 0;  //suspisiouce - will be changed
 777
 778     if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)
 779         // strict cut on overlaped cluster
 780         return 0;  //suspisiouce - will be changed
 781
 782     if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 )
 783          && t.fCurrentCluster->GetType()<0)
 784       return 0;
 785
 786     //    t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
 787     UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
 788
 789   } else {
 790     if (y > ymax) {
 791       t.fRelativeSector= (t.fRelativeSector+1) % fN;
 792       if (!t.Rotate(fSectors->GetAlpha()))
 793         return 0;
 794     } else if (y <-ymax) {
 795       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
 796       if (!t.Rotate(-fSectors->GetAlpha()))
 797         return 0;
 798     }
 799   }
 800   return 1;
 801 }
 802
 803
 804 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex,  Int_t nr) {
 805   //-----------------------------------------------------------------
 806   // This function tries to find a track prolongation to next pad row
 807   //-----------------------------------------------------------------
 808   t.fCurrentCluster  = 0;
 809   t.fCurrentClusterIndex1 = 0;
 810   t.fCurrentClusterIndex2 = 0;
 811
 812   Double_t xt=t.GetX();
 813   Int_t row = fSectors->GetRowNumber(xt)-1;
 814   if (row < nr) return 1; // don't prolongate if not information until now -
 815   Double_t x=fSectors->GetX(nr);
 816   //  if (t.fStopped) return 0;
 817   //  if (t.GetRadius()>x+10 ) return 0;
 818   if (!t.PropagateTo(x)){
 819     t.fStopped =kTRUE;
 820     return 0;
 821   }
 822   // update current
 823   t.fCurrentSigmaY = GetSigmaY(&t);
 824   t.fCurrentSigmaZ = GetSigmaZ(&t);
 825
 826   AliTPCclusterMI *cl=0;
 827   UInt_t index=0;
 828   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
 829   //
 830   Double_t  y=t.GetY(), z=t.GetZ();
 831   Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
 832   Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
 833   //
 834
 835   Float_t maxdistance = 1000000;
 836   if (krow) {
 837     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
 838       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
 839       if (c->GetZ() > z+roadz) break;
 840       if (TMath::Abs(c->GetY()-y)>roady) continue;
 841
 842       //krow.UpdateClusterTrack(i,trindex,&t);
 843
 844       Float_t dy2 = (c->GetY()- t.GetY());
 845       dy2*=dy2;
 846       Float_t dz2 = (c->GetZ()- t.GetZ());
 847       dz2*=dz2;
 848       //
 849       Float_t distance = dy2+dz2;
 850       //
 851       if (distance > maxdistance) continue;
 852       maxdistance = distance;
 853       cl=c;
 854       index=i;
 855     }
 856   }
 857   t.fCurrentCluster  = cl;
 858   t.fCurrentClusterIndex1 = krow.GetIndex(index);
 859   t.fCurrentClusterIndex2 = index;
 860   return 1;
 861 }
 862
 863
 864 Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
 865   //-----------------------------------------------------------------
 866   // This function tries to find a track prolongation to next pad row
 867   //-----------------------------------------------------------------
 868   AliTPCseed & t  = *((AliTPCseed*)(fSeeds->At(trindex)));
 869   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
 870   //  Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
 871   Double_t y=t.GetY();
 872   Double_t ymax=fSectors->GetMaxY(nr);
 873
 874   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
 875     t.fInDead = kTRUE;
 876     Int_t row = nr;
 877     if (fSectors==fOuterSec) row += fParam->GetNRowLow();
 878     t.fClusterIndex[row] = -1;
 879     return 0;
 880   }
 881   else
 882     {
 883       if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
 884       else
 885         return 0;
 886     }
 887
 888   if (t.fCurrentCluster) {
 889     //    Float_t l=fSectors->GetPadPitchLength();
 890     //    AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
 891
 892     Double_t sy2=ErrY2(&t,t.fCurrentCluster);
 893     Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
 894
 895
 896     Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
 897     Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
 898
 899     Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
 900     Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
 901
 902     Double_t rdistance  = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
 903
 904
 905     //    printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
 906     if ( (rdistancey>1) || (rdistancez>1)) return 0;
 907     if (rdistance>4) return 0;
 908
 909     if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)
 910         return 0;  //suspisiouce - will be changed
 911
 912     if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)
 913         // strict cut on overlaped cluster
 914         return 0;  //suspisiouce - will be changed
 915
 916     if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 )
 917          && t.fCurrentCluster->GetType()<0)
 918       return 0;
 919
 920     //    t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
 921     UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
 922
 923   } else {
 924     if (y > ymax) {
 925       t.fRelativeSector= (t.fRelativeSector+1) % fN;
 926       if (!t.Rotate(fSectors->GetAlpha()))
 927         return 0;
 928     } else if (y <-ymax) {
 929       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
 930       if (!t.Rotate(-fSectors->GetAlpha()))
 931         return 0;
 932     }
 933   }
 934   return 1;
 935 }
 936
 937
 938
 939
 940 /*
 941 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
 942 {
 943   //-----------------------------------------------------------------
 944   // fast prolongation mathod -
 945   // don't update track only after step clusters
 946   //-----------------------------------------------------------------
 947   Double_t xt=t.GetX();
 948   //
 949   Double_t alpha=t.GetAlpha();
 950   alpha =- fSectors->GetAlphaShift();
 951   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 952   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 953   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
 954   Int_t row0 = fSectors->GetRowNumber(xt);
 955   Double_t x    = fSectors->GetX(row0);
 956   Double_t ymax = fSectors->GetMaxY(row0);
 957   //
 958   Double_t sy2=ErrY2(&t)*2;
 959   Double_t sz2=ErrZ2(&t)*2;
 960   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
 961   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
 962   Float_t maxdistance = roady*roady + roadz*roadz;
 963   t.fCurrentSigmaY = GetSigmaY(&t);
 964   t.fCurrentSigmaZ = GetSigmaZ(&t);
 965   //
 966   Int_t nclusters = 0;
 967   Double_t y;
 968   Double_t z;
 969   Double_t yy[200];   //track prolongation
 970   Double_t zz[200];
 971   Double_t cy[200];  // founded cluster position
 972   Double_t cz[200];
 973   Double_t sy[200];  // founded cluster error
 974   Double_t sz[200];
 975   Bool_t   hitted[200]; // indication of cluster presence
 976   //
 977
 978   //
 979   for (Int_t drow = step; drow>=0; drow--) {
 980     Int_t row = row0-drow;
 981     if (row<0) break;
 982     Double_t x    = fSectors->GetX(row);
 983     Double_t ymax = fSectors->GetMaxY(row);
 984     t.GetProlongation(x,y,z);
 985     yy[drow] =y;
 986     zz[drow] =z;
 987     const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
 988     if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
 989       t.fInDead = kTRUE;
 990       break;
 991     }
 992     else
 993       {
 994         t.fNFoundable++;
 995       }
 996
 997     //find nearest  cluster
 998     AliTPCclusterMI *cl= 0;
 999     if (krow) {
1000       for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
1001         AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
1002         if (c->GetZ() > z+roadz) break;
1003         if ( (c->GetY()-y) >  roady ) continue;
1004         Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
1005         if (maxdistance>distance) {
1006           maxdistance = distance;
1007           cl=c;
1008           //      index=krow.GetIndex(i);
1009         }
1010       }
1011     }  // end of seearch
1012     //update cluster information
1013     if (cl){
1014       cy[drow] = cl->GetY();
1015       cz[drow] = cl->GetZ();
1016       sy[drow] = ErrY2(&t,cl);
1017       sz[drow] = ErrZ2(&t,cl);
1018       hitted[drow] = kTRUE;
1019       nclusters++;
1020     }
1021     else
1022       hitted[drow] = kFALSE;
1023   }
1024   //if we have information - update track
1025   if (nclusters>0){
1026     Float_t sumyw0   = 0;
1027     Float_t sumzw0   = 0;
1028     Float_t sumyw   = 0;
1029     Float_t sumzw   = 0;
1030     Float_t sumrow  = 0;
1031     for (Int_t i=0;i<step;i++)
1032       if (hitted[i]){
1033         sumyw0+= 1/sy[i];
1034         sumzw0+= 1/sz[i];
1035         //
1036         sumyw+= (cy[i]-yy[i])/sy[i];
1037         sumzw+= (cz[i]-zz[i])/sz[i];
1038         sumrow+= i;
1039       }
1040     Float_t dy = sumyw/sumyw0;
1041     Float_t dz = sumzw/sumzw0;
1042     Float_t mrow = sumrow/nclusters+row0;
1043     Float_t x = fSectors->GetX(mrow);
1044     t.PropagateTo(x);
1045     AliTPCclusterMI cvirtual;
1046     cvirtual.SetZ(dz+t.GetZ());
1047     cvirtual.SetY(dy+t.GetY());
1048     t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
1049     t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
1050     Float_t chi2 = t.GetPredictedChi2(&cvirtual);
1051     t.Update(&cvirtual,chi2,0);
1052     Int_t ncl = t.GetNumberOfClusters();
1053     ncl = ncl-1+nclusters;
1054     t.SetN(ncl);
1055   }
1056   return  1;
1057 }
1058 */
1059
1060 //_____________________________________________________________________________
1061 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
1062   //-----------------------------------------------------------------
1063   // This function tries to find a track prolongation.
1064   //-----------------------------------------------------------------
1065   Double_t xt=t.GetX();
1066   //
1067   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1068   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1069   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1070   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
1071
1072   for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
1073
1074     if (FollowToNext(t,nr)==0) {
1075     }
1076   }
1077   return 1;
1078 }
1079
1080
1081 //_____________________________________________________________________________
1082 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1083   //-----------------------------------------------------------------
1084   // This function tries to find a track prolongation.
1085   //-----------------------------------------------------------------
1086   Double_t xt=t.GetX();
1087   //
1088   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1089   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1090   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1091   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
1092
1093   for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
1094     if (t.GetSnp()<0.8)
1095       FollowToNext(t,nr);
1096   }
1097   return 1;
1098 }
1099
1100
1101
1102
1103
1104 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1105 {
1106   //
1107   //
1108   sum1=0;
1109   sum2=0;
1110   Int_t sum=0;
1111   //
1112   Float_t dz2 =(s1->GetZ() - s2->GetZ());
1113   dz2*=dz2;
1114   /*
1115     Float_t x = s1->GetX();
1116     Float_t x2 = s2->GetX();
1117
1118     Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
1119   */
1120   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1121   //if (TMath::Abs(dy2)>2*ymax-3)
1122   //  dy2-=2*ymax;
1123   dy2*=dy2;
1124   Float_t distance = TMath::Sqrt(dz2+dy2);
1125   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
1126
1127   Int_t offset =0;
1128   if (fSectors==fOuterSec) offset = fParam->GetNRowLow();
1129   Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
1130   Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
1131   lastpoint +=offset;
1132   firstpoint+=offset;
1133   if (lastpoint>160)
1134     lastpoint =160;
1135   if (firstpoint<0)
1136     firstpoint = 0;
1137   if (firstpoint<lastpoint-15) {
1138     firstpoint =0;
1139     lastpoint  =160;
1140   }
1141
1142
1143   for (Int_t i=firstpoint;i<lastpoint;i++){
1144     if (s1->fClusterIndex[i]>0) sum1++;
1145     if (s2->fClusterIndex[i]>0) sum2++;
1146     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
1147       sum++;
1148     }
1149   }
1150
1151   Float_t summin = TMath::Min(sum1+1,sum2+1);
1152   Float_t ratio = (sum+1)/Float_t(summin);
1153   return ratio;
1154 }
1155
1156 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1157 {
1158   //
1159   //
1160   if (s1->fSector!=s2->fSector) return;
1161   //
1162   Float_t dz2 =(s1->GetZ() - s2->GetZ());
1163   dz2*=dz2;
1164   Float_t dy2 =(s1->GetY() - s2->GetY());
1165
1166   dy2*=dy2;
1167   Float_t distance = TMath::Sqrt(dz2+dy2);
1168   if (distance>15.) return ; // if there are far away  - not overlap - to reduce combinatorics
1169   //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
1170   //  TClonesArray &pointarray1 = *(s1->fPoints);
1171   //TClonesArray &pointarray2 = *(s2->fPoints);
1172   //
1173   for (Int_t i=0;i<160;i++){
1174     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
1175       //  AliTPCTrackPoint *p1  = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
1176       //AliTPCTrackPoint *p2  = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i));
1177       AliTPCTrackPoint *p1  = s1->GetTrackPoint(i);
1178       AliTPCTrackPoint *p2  = s2->GetTrackPoint(i);;
1179       p1->fIsShared = kTRUE;
1180       p2->fIsShared = kTRUE;
1181     }
1182   }
1183 }
1184
1185
1186
1187
1188 void  AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared){
1189
1190
1191
1192   //
1193   // remove overlap - used removal factor - removal index stored in the track
1194   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1195     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1196     if (pt) RotateToLocal(pt);
1197   }
1198   arr->Sort();  // sorting according z
1199   arr->Expand(arr->GetEntries());
1200   Int_t nseed=arr->GetEntriesFast();
1201   //  printf("seeds \t%p \t%d\n",arr, nseed);
1202   //  arr->Expand(arr->GetEntries());  //remove 0 pointers
1203   nseed = arr->GetEntriesFast();
1204   Int_t removed = 0;
1205   for (Int_t i=0; i<nseed; i++) {
1206     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1207     if (!pt) {
1208       continue;
1209     }
1210     if (!(pt->IsActive())) continue;
1211     for (Int_t j=i+1; j<nseed; j++){
1212       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1213       //
1214       if (!pt2) continue;
1215       if (!(pt2->IsActive())) continue;
1216       if (TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
1217       if (TMath::Abs(pt2->GetZ()-pt->GetZ())<4){
1218         Int_t sum1,sum2;
1219         Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
1220         //if (sum1==0) {
1221         //  pt->Desactivate(removalindex); // arr->RemoveAt(i);
1222         //  break;
1223         //}
1224         if (ratio>factor){
1225           //      if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex);  // arr->RemoveAt(j);
1226           Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
1227           Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
1228           removed++;
1229           if (TMath::Abs(ratio3)>0.025){  // if much more points
1230             if (sum1>sum2) pt2->Desactivate(removalindex);
1231             else {
1232               pt->Desactivate(removalindex); // arr->RemoveAt(i);
1233               break;
1234             }
1235           }
1236           else{  //decide on mean chi2
1237             if (ratio2<1)
1238               pt2->Desactivate(removalindex);
1239             else {
1240               pt->Desactivate(removalindex); // arr->RemoveAt(i);
1241               break;
1242             }
1243           }
1244
1245         }  // if suspicious ratio
1246       }
1247       else
1248         break;
1249     }
1250   }
1251   //  printf("removed\t%d\n",removed);
1252   Int_t good =0;
1253   for (Int_t i=0; i<nseed; i++) {
1254     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1255     if (!pt) break;
1256     if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
1257       //desactivate tracks with small number of points
1258       //      printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
1259       pt->Desactivate(10);  //desactivate  - small muber of points
1260     }
1261     if (!(pt->IsActive())) continue;
1262     good++;
1263   }
1264
1265
1266   if (shared)
1267     for (Int_t i=0; i<nseed; i++) {
1268       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1269       if (!pt) continue;
1270       if (!(pt->IsActive())) continue;
1271       for (Int_t j=i+1; j<nseed; j++){
1272         AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1273         if ((pt2) && pt2->IsActive()) {
1274           if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
1275           SignShared(pt,pt2);
1276         }
1277       }
1278     }
1279   fNtracks = good;
1280   printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
1281
1282
1283 }
1284
1285 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor, Int_t removalindex)
1286 {
1287
1288   //Loop over all tracks and remove "overlaps"
1289   //
1290   //
1291   Int_t nseed = arr->GetEntriesFast();
1292   Int_t good =0;
1293   for (Int_t i=0; i<nseed; i++) {
1294     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1295     if (!pt) {
1296       continue;
1297     }
1298     if (!(pt->IsActive())) continue;
1299     Int_t noc=pt->GetNumberOfClusters();
1300     Int_t shared =0;
1301     for (Int_t i=0; i<noc; i++) {
1302       Int_t index=pt->GetClusterIndex(i);
1303       AliTPCclusterMI *c=(AliTPCclusterMI*)GetClusterMI(index);
1304       if (!c) continue;
1305       if (c->IsUsed()) shared++;
1306     }
1307     if ((Float_t(shared)/Float_t(noc))>factor)
1308       pt->Desactivate(removalindex);
1309     else{
1310       good++;
1311       for (Int_t i=0; i<noc; i++) {
1312         Int_t index=pt->GetClusterIndex(i);
1313         AliTPCclusterMI *c=(AliTPCclusterMI*)GetClusterMI(index);
1314         if (!c) continue;
1315         c->Use();
1316       }
1317     }
1318   }
1319   fNtracks = good;
1320   printf("\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
1321
1322 }
1323
1324
1325 void AliTPCtrackerMI::MakeSeedsAll()
1326 {
1327   if (fSeeds == 0) fSeeds = new TObjArray;
1328   TObjArray * arr;
1329   for (Int_t sec=0;sec<fkNOS;sec+=3){
1330      arr = MakeSeedsSectors(sec,sec+3);
1331      Int_t nseed = arr->GetEntriesFast();
1332      for (Int_t i=0;i<nseed;i++)
1333        fSeeds->AddLast(arr->RemoveAt(i));
1334   }
1335   //  fSeeds = MakeSeedsSectors(0,fkNOS);
1336 }
1337
1338 TObjArray *  AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
1339 {
1340   //
1341   // loop over all  sectors and make seed
1342   //find track seeds
1343   TStopwatch timer;
1344   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1345   Int_t nrows=nlow+nup;
1346   Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
1347   //  if (fSeeds==0) fSeeds = new TObjArray;
1348   TObjArray * arr = new TObjArray;
1349
1350   for (Int_t sec=sec1; sec<sec2;sec++){
1351     MakeSeeds(arr, sec, nup-1, nup-1-gap);
1352     MakeSeeds(arr, sec, nup-2-shift, nup-2-shift-gap);
1353   }
1354   gap = Int_t(0.2* nrows);
1355   for (Int_t sec=sec1; sec<sec2;sec++){
1356     //find secondaries
1357     //MakeSeeds2(arr, sec, nup-1, nup-1-gap);
1358     MakeSeeds2(arr, sec, nup-1-shift, nup-1-shift-gap);
1359     //MakeSeeds2(arr, sec, nup-1-2*shift, nup-1-2*shift-gap);
1360     MakeSeeds2(arr, sec, nup-1-3*shift, nup-1-3*shift-gap);
1361     //MakeSeeds2(arr, sec, nup-1-4*shift, nup-1-4*shift-gap);
1362     MakeSeeds2(arr, sec, nup-1-5*shift, nup-1-5*shift-gap);
1363     MakeSeeds2(arr, sec, gap, 1);
1364   }
1365
1366   Int_t nseed=arr->GetEntriesFast();
1367   Int_t i;
1368
1369     gap=Int_t(0.3*nrows);
1370     // continue seeds
1371     for (i=0; i<nseed; i++) {
1372     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1373     if (!pt) continue;
1374     if (FollowProlongation(t,nup-gap)) {
1375     pt->fIsSeeding =kFALSE;
1376     continue;
1377     }
1378     delete arr->RemoveAt(i);
1379     }
1380
1381
1382   //
1383   //remove seeds which overlaps
1384   RemoveOverlap(arr,0.4,1);
1385   //delete seeds - which were sign
1386   nseed=arr->GetEntriesFast();
1387   for (i=0; i<nseed; i++) {
1388     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1389     //, &t=*pt;
1390     if (!pt) continue;
1391     if ((pt->IsActive())  &&  pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
1392       //pt->Reset();
1393       //FollowBackProlongation(*pt,nup-1);
1394       //if ( pt->GetNumberOfClusters() < pt->fNFoundable*0.5 || pt->GetNumberOfClusters()<10 )
1395       //delete arr->RemoveAt(i);
1396       //else
1397       //        pt->Reset();
1398       continue;
1399     }
1400     delete arr->RemoveAt(i);
1401   }
1402   //RemoveOverlap(arr,0.6,1);
1403   return arr;
1404 }
1405
1406
1407
1408 //_____________________________________________________________________________
1409 void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1410   //-----------------------------------------------------------------
1411   // This function creates track seeds.
1412   //-----------------------------------------------------------------
1413   //  if (fSeeds==0) fSeeds=new TObjArray(15000);
1414
1415   Double_t x[5], c[15];
1416
1417   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1418   Double_t cs=cos(alpha), sn=sin(alpha);
1419
1420   Double_t x1 =fOuterSec->GetX(i1);
1421   Double_t xx2=fOuterSec->GetX(i2);
1422   //
1423   //  for (Int_t ns=0; ns<fkNOS; ns++)
1424   Int_t ns =sec;
1425     {
1426     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1427     Int_t nm=fOuterSec[ns][i2];
1428     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
1429     const AliTPCRow& kr1=fOuterSec[ns][i1];
1430     AliTPCRow&  kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1431     AliTPCRow&  kr22 = fOuterSec[(ns)%fkNOS][i2];
1432     AliTPCRow&  kr23 = fOuterSec[(ns+1)%fkNOS][i2];
1433
1434     for (Int_t is=0; is < kr1; is++) {
1435       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1436       Double_t x3=GetX(), y3=GetY(), z3=GetZ();
1437
1438       Float_t anglez = (z1-z3)/(x1-x3);
1439       Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
1440
1441       for (Int_t js=0; js < nl+nm+nu; js++) {
1442         const AliTPCclusterMI *kcl;
1443         Double_t x2,   y2,   z2;
1444         if (js<nl) {
1445           if (js==0) {
1446             js = kr21.Find(extraz-15.);
1447             if (js>=nl) continue;
1448           }
1449           kcl=kr21[js];
1450           z2=kcl->GetZ();
1451           if ((extraz-z2)>10) continue;
1452           if ((extraz-z2)<-10) {
1453             js = nl-1;
1454             continue;
1455           }
1456           y2=kcl->GetY();
1457           x2= xx2*cs+y2*sn;
1458           y2=-xx2*sn+y2*cs;
1459         } else
1460           if (js<nl+nm) {
1461             if (js==nl) {
1462               js = nl+kr22.Find(extraz-15.);
1463               if (js>=nl+nm) continue;
1464             }
1465             kcl=kr22[js-nl];
1466             z2=kcl->GetZ();
1467             if ((extraz-z2)>10) continue;
1468             if ((extraz-z2)<-10) {
1469               js = nl+nm-1;
1470               continue;
1471             }
1472             x2=xx2; y2=kcl->GetY();
1473           } else {
1474             //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
1475             if (js==nl+nm) {
1476               js = nl+nm+kr23.Find(extraz-15.);
1477               if (js>=nl+nm+nu) break;
1478             }
1479             kcl=kr23[js-nl-nm];
1480             z2=kcl->GetZ();
1481             if ((extraz-z2)>10) continue;
1482             if ((extraz-z2)<-10) {
1483               break;
1484             }
1485             y2=kcl->GetY();
1486             x2=xx2*cs-y2*sn;
1487             y2=xx2*sn+y2*cs;
1488           }
1489
1490         Double_t zz=z1 - anglez*(x1-x2);
1491         if (TMath::Abs(zz-z2)>10.) continue;
1492
1493         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1494         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
1495
1496         x[0]=y1;
1497         x[1]=z1;
1498         x[4]=f1(x1,y1,x2,y2,x3,y3);
1499         if (TMath::Abs(x[4]) >= 0.0066) continue;
1500         x[2]=f2(x1,y1,x2,y2,x3,y3);
1501         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1502         x[3]=f3(x1,y1,x2,y2,z1,z2);
1503         if (TMath::Abs(x[3]) > 1.2) continue;
1504         Double_t a=asin(x[2]);
1505         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1506         if (TMath::Abs(zv-z3)>10.) continue;
1507
1508         Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
1509         Double_t sy2=kcl->GetSigmaY2()*2,     sz2=kcl->GetSigmaZ2()*4;
1510         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1511         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1512         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1513
1514         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1515         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1516         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1517         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1518         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1519         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1520         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1521         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1522         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1523         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1524
1525         c[0]=sy1;
1526         c[1]=0.;       c[2]=sz1;
1527         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1528         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1529                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1530         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1531         c[13]=f30*sy1*f40+f32*sy2*f42;
1532         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1533
1534         UInt_t index=kr1.GetIndex(is);
1535         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
1536         track->fIsSeeding = kTRUE;
1537         //Int_t rc=FollowProlongation(*track, i2);
1538         Int_t delta4 = Int_t((i2-i1)/4.);
1539
1540         FollowProlongation(*track, i1-delta4);
1541         if (track->GetNumberOfClusters() < track->fNFoundable/2.) {
1542           delete track;
1543           continue;
1544         }
1545         FollowProlongation(*track, i1-2*delta4);
1546         if (track->GetNumberOfClusters() < track->fNFoundable/2.) {
1547           delete track;
1548           continue;
1549         }
1550         FollowProlongation(*track, i1-3*delta4);
1551         if (track->GetNumberOfClusters() < track->fNFoundable/2.) {
1552           delete track;
1553           continue;
1554         }
1555         FollowProlongation(*track, i2);
1556         //Int_t rc = 1;
1557
1558         track->fLastPoint = i1;  // first cluster in track position
1559         if (track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1560         else arr->AddLast(track);
1561       }
1562     }
1563   }
1564 }
1565
1566
1567 //_____________________________________________________________________________
1568 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1569   //-----------------------------------------------------------------
1570   // This function creates track seeds - without vertex constraint
1571   //-----------------------------------------------------------------
1572
1573   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1574   //  Double_t cs=cos(alpha), sn=sin(alpha);
1575   Int_t row0 = (i1+i2)/2;
1576   Int_t drow = (i1-i2)/2;
1577   const AliTPCRow& kr0=fSectors[sec][row0];
1578   const AliTPCRow& krm=fSectors[sec][row0-1];
1579   const AliTPCRow& krp=fSectors[sec][row0+1];
1580   AliTPCRow * kr=0;
1581
1582   AliTPCpolyTrack polytrack;
1583   Int_t nclusters=fSectors[sec][row0];
1584
1585   for (Int_t is=0; is < nclusters; is++) {
1586     const AliTPCclusterMI * cl= kr0[is];
1587     Double_t x = kr0.GetX();
1588
1589     // Initialization of the polytrack
1590     polytrack.Reset();
1591
1592     Double_t y0= cl->GetY();
1593     Double_t z0= cl->GetZ();
1594     polytrack.AddPoint(x,y0,z0);
1595     Float_t roady = 5*TMath::Sqrt(cl->GetSigmaY2()+0.2);
1596     Float_t roadz = 5*TMath::Sqrt(cl->GetSigmaZ2()+0.2);
1597     //
1598     x = krm.GetX();
1599     cl = krm.FindNearest(y0,z0,roady,roadz);
1600     if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),roady,roadz);
1601     //
1602     x = krp.GetX();
1603     cl = krp.FindNearest(y0,z0,roady,roadz);
1604     if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),cl->GetSigmaY2()+0.05,cl->GetSigmaZ2()+0.05);
1605     //
1606     polytrack.UpdateParameters();
1607     // follow polytrack
1608     roadz = 0.6;
1609     roady = 0.6;
1610     //
1611     Double_t yn,zn;
1612     Int_t nfoundable = polytrack.GetN();
1613     Int_t nfound     = nfoundable;
1614     for (Int_t ddrow = 2; ddrow<drow;ddrow++){
1615       for (Int_t delta = -1;delta<=1;delta+=2){
1616         Int_t row = row0+ddrow*delta;
1617         kr = &(fSectors[sec][row]);
1618         Double_t xn = kr->GetX();
1619         Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone;
1620         polytrack.GetFitPoint(xn,yn,zn);
1621         if (TMath::Abs(yn)>ymax) continue;
1622         nfoundable++;
1623         AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
1624         if (cln) {
1625           polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),cln->GetSigmaY2()+0.05,cln->GetSigmaZ2()+0.05);
1626           nfound++;
1627         }
1628       }
1629       polytrack.UpdateParameters();
1630       if (nfound<0.45*nfoundable) break;
1631     }
1632     if ((nfound>0.5*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
1633       // add polytrack candidate
1634       Double_t x[5], c[15];
1635       Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
1636       polytrack.GetBoundaries(x3,x1);
1637       x2 = (x1+x3)/2.;
1638       polytrack.GetFitPoint(x1,y1,z1);
1639       polytrack.GetFitPoint(x2,y2,z2);
1640       polytrack.GetFitPoint(x3,y3,z3);
1641       //
1642       //is track pointing to the vertex ?
1643       Double_t x0,y0,z0;
1644       x0=0;
1645       polytrack.GetFitPoint(x0,y0,z0);
1646       if ( (TMath::Abs(z0-GetZ())<10) && (TMath::Abs(y0-GetY())<5)){ //if yes apply vertex constraint
1647         //      x3 = 0;
1648         //y3 = GetY();
1649         //z3 = GetZ();
1650       }
1651
1652       x[0]=y1;
1653       x[1]=z1;
1654       x[4]=f1(x1,y1,x2,y2,x3,y3);
1655       if (TMath::Abs(x[4]) >= 0.05) continue;  //MI change
1656       x[2]=f2(x1,y1,x2,y2,x3,y3);
1657       //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1658       x[3]=f3(x1,y1,x2,y2,z1,z2);
1659       if (TMath::Abs(x[3]) > 1.2) continue;
1660       if (TMath::Abs(x[2]) > 0.99) continue;
1661       //      Double_t a=asin(x[2]);
1662
1663
1664       Double_t sy=1.5, sz=1.5;
1665       Double_t sy1=1.5, sz1=1.5;
1666       Double_t sy2=1.3, sz2=1.3;
1667       Double_t sy3=1.5;
1668       //sz3=1.5;
1669       //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1670       //      Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1671       //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1672
1673       Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1674       Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1675       Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1676       Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1677       Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1678       Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1679       Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1680       Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1681       Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1682       Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1683
1684       c[0]=sy1;
1685       c[1]=0.;       c[2]=sz1;
1686       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1687       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1688       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1689       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1690       c[13]=f30*sy1*f40+f32*sy2*f42;
1691       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1692
1693       UInt_t index=0;
1694       //kr0.GetIndex(is);
1695       AliTPCseed *track=new AliTPCseed(index, x, c, x1, sec*alpha+shift);
1696       track->fStopped =kFALSE;
1697       track->fIsSeeding = kTRUE;
1698       Int_t rc=FollowProlongation(*track, i2);
1699         track->fLastPoint = i1;  // first cluster in track position
1700         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1701         else arr->AddLast(track);
1702     }
1703   }
1704
1705
1706
1707 }
1708
1709
1710
1711
1712
1713 //_____________________________________________________________________________
1714 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
1715   //-----------------------------------------------------------------
1716   // This function reades track seeds.
1717   //-----------------------------------------------------------------
1718   TDirectory *savedir=gDirectory;
1719
1720   TFile *in=(TFile*)inp;
1721   if (!in->IsOpen()) {
1722      cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
1723      return 1;
1724   }
1725
1726   in->cd();
1727   TTree *seedTree=(TTree*)in->Get("Seeds");
1728   if (!seedTree) {
1729      cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
1730      cerr<<"can't get a tree with track seeds !\n";
1731      return 2;
1732   }
1733   AliTPCtrack *seed=new AliTPCtrack;
1734   seedTree->SetBranchAddress("tracks",&seed);
1735
1736   if (fSeeds==0) fSeeds=new TObjArray(15000);
1737
1738   Int_t n=(Int_t)seedTree->GetEntries();
1739   for (Int_t i=0; i<n; i++) {
1740      seedTree->GetEvent(i);
1741      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
1742   }
1743
1744   delete seed;
1745   delete seedTree;
1746   savedir->cd();
1747   return 0;
1748 }
1749
1750 //_____________________________________________________________________________
1751 Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
1752   //-----------------------------------------------------------------
1753   // This is a track finder.
1754   //-----------------------------------------------------------------
1755   TDirectory *savedir=gDirectory;
1756
1757   if (inp) {
1758      TFile *in=(TFile*)inp;
1759      if (!in->IsOpen()) {
1760         cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
1761         return 1;
1762      }
1763   }
1764
1765   if (!out->IsOpen()) {
1766      cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
1767      return 2;
1768   }
1769
1770   out->cd();
1771
1772   char   tname[100];
1773   sprintf(tname,"TreeT_TPC_%d",fEventN);
1774   TTree tracktree(tname,"Tree with TPC tracks");
1775   TTree seedtree("Seeds","Seeds");
1776   AliTPCtrack *iotrack=0;
1777   AliTPCseed  *ioseed=0;
1778   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
1779   TStopwatch timer;
1780
1781   printf("Loading clusters \n");
1782   LoadClusters();
1783   printf("Time for loading clusters: \t");timer.Print();timer.Start();
1784
1785   printf("Loading outer sectors\n");
1786   LoadOuterSectors();
1787   printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
1788
1789   printf("Loading inner sectors\n");
1790   LoadInnerSectors();
1791   printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
1792   fSectors = fOuterSec;
1793   fN=fkNOS;
1794
1795
1796   //find track seeds
1797   MakeSeedsAll();
1798   printf("Time for seeding: \t"); timer.Print();timer.Start();
1799   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1800   Int_t nrows=nlow+nup;
1801
1802   Int_t gap=Int_t(0.3*nrows);
1803   Int_t i;
1804   //RemoveOverlap(fSeeds,0.6,2);
1805   Int_t nseed=fSeeds->GetEntriesFast();
1806   // outer sectors parallel tracking
1807   ParallelTracking(fSectors->GetNRows()-gap-1,0);
1808   printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
1809
1810   RemoveOverlap(fSeeds, 0.4,3);
1811   printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
1812   //parallel tracking
1813   fSectors = fInnerSec;
1814   fN=fkNIS;
1815
1816   ParallelTracking(fSectors->GetNRows()-1,0);
1817   /*
1818     ParallelTracking(fSectors->GetNRows()-1,2*fSectors->GetNRows()/3);
1819     RemoveOverlap(fSeeds,0.4,5,kTRUE);
1820     ParallelTracking(2*fSectors->GetNRows()/3-1,fSectors->GetNRows()/3);
1821     RemoveOverlap(fSeeds,0.4,5,kTRUE);
1822     ParallelTracking(fSectors->GetNRows()/3-1,0);
1823   */
1824   printf("Number of tracks after  inner tracking  %d\n",fNtracks);
1825   printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
1826   //
1827   for (Int_t i=0;i<fSeeds->GetEntriesFast();i++){
1828     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
1829     if (!pt) continue;
1830     if (!pt->IsActive()) continue;
1831     pt->PropagateTo(90.);
1832   }
1833   RemoveOverlap(fSeeds,0.4,5,kTRUE);  // remove overlap -  shared points signed
1834   RemoveUsed(fSeeds,0.4,6);
1835   printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
1836   //
1837   //
1838
1839   ioseed  = (AliTPCseed*)(fSeeds->UncheckedAt(0));
1840   AliTPCseed * vseed = new AliTPCseed;
1841   vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
1842   vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
1843   vseed->fPoints->ExpandCreateFast(2);
1844
1845   //TBranch * seedbranch =
1846   seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
1847   //delete vseed;
1848   nseed=fSeeds->GetEntriesFast();
1849
1850   Int_t found = 0;
1851   for (i=0; i<nseed; i++) {
1852     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1853     if (!pt) continue;
1854     Int_t nc=t.GetNumberOfClusters();
1855     if (nc<20) continue;
1856     t.CookdEdx(0.02,0.6);
1857     CookLabel(pt,0.1); //For comparison only
1858     //   if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
1859     if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (t.fNFoundable>Int_t(0.3*nrows)))){
1860       iotrack=pt;
1861       tracktree.Fill();
1862      cerr<<found++<<'\r';
1863     }
1864     /*
1865       pt->RebuildSeed();
1866       seedbranch->SetAddress(&pt);
1867
1868       seedtree.Fill();
1869       for (Int_t j=0;j<160;j++){
1870       delete pt->fPoints->RemoveAt(j);
1871       }
1872       delete pt->fPoints;
1873       pt->fPoints =0;
1874     */
1875     delete fSeeds->RemoveAt(i);
1876   }
1877   //  fNTracks = found;
1878   printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
1879
1880   UnloadClusters();
1881   printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1882
1883   tracktree.Write();
1884   seedtree.Write();
1885   cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
1886
1887   savedir->cd();
1888
1889   return 0;
1890 }
1891
1892
1893 void  AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
1894 {
1895   //
1896   // try to track in parralel
1897
1898   Int_t nseed=fSeeds->GetEntriesFast();
1899   //prepare seeds for tracking
1900   for (Int_t i=0; i<nseed; i++) {
1901     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1902     if (!pt) continue;
1903     if (!t.IsActive()) continue;
1904     // follow prolongation to the first layer
1905     if ( (fSectors ==fInnerSec) || (t.fFirstPoint>rfirst+1))
1906       FollowProlongation(t, rfirst+1);
1907   }
1908
1909
1910   //
1911   for (Int_t nr=rfirst; nr>=rlast; nr--){
1912     // make indexes with the cluster tracks for given
1913     //    for (Int_t i = 0;i<fN;i++)
1914     //  fSectors[i][nr].MakeClusterTracks();
1915
1916     // find nearest cluster
1917     for (Int_t i=0; i<nseed; i++) {
1918       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1919       if (!pt) continue;
1920       if (!pt->IsActive()) continue;
1921       if ( (fSectors ==fOuterSec) && pt->fFirstPoint<nr) continue;
1922       if (pt->fRelativeSector>17) {
1923         continue;
1924       }
1925       UpdateClusters(t,i,nr);
1926     }
1927     // prolonagate to the nearest cluster - if founded
1928     for (Int_t i=0; i<nseed; i++) {
1929       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
1930       if (!pt) continue;
1931       if (!pt->IsActive()) continue;
1932       if ((fSectors ==fOuterSec) &&pt->fFirstPoint<nr) continue;
1933       if (pt->fRelativeSector>17) {
1934         continue;
1935       }
1936       FollowToNextCluster(i,nr);
1937     }
1938     //    for (Int_t i= 0;i<fN;i++)
1939     //  fSectors[i][nr].ClearClusterTracks();
1940   }
1941
1942 }
1943
1944 Float_t  AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
1945 {
1946   //
1947   //
1948   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
1949   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1950   Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
1951   Float_t angular  = seed->GetSnp();
1952   angular = angular*angular/(1-angular*angular);
1953   //  angular*=angular;
1954   //angular  = TMath::Sqrt(angular/(1-angular));
1955   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
1956   return res;
1957 }
1958 Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
1959 {
1960   //
1961   //
1962   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
1963   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1964   Float_t sres = fParam->GetZSigma();
1965   Float_t angular  = seed->GetTgl();
1966   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
1967   return res;
1968 }
1969
1970
1971
1972 //_________________________________________________________________________
1973 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1974   //--------------------------------------------------------------------
1975   //       Return pointer to a given cluster
1976   //--------------------------------------------------------------------
1977   Int_t sec=(index&0xff000000)>>24;
1978   Int_t row=(index&0x00ff0000)>>16;
1979   Int_t ncl=(index&0x0000ffff)>>00;
1980
1981   AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
1982   if (!clrow) return 0;
1983   return (AliTPCclusterMI*)(*clrow)[ncl];
1984 }
1985
1986 //__________________________________________________________________________
1987 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1988   //--------------------------------------------------------------------
1989   //This function "cooks" a track label. If label<0, this track is fake.
1990   //--------------------------------------------------------------------
1991   Int_t noc=t->GetNumberOfClusters();
1992   Int_t *lb=new Int_t[noc];
1993   Int_t *mx=new Int_t[noc];
1994   AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
1995
1996   Int_t i;
1997   for (i=0; i<noc; i++) {
1998      lb[i]=mx[i]=0;
1999      Int_t index=t->GetClusterIndex(i);
2000      clusters[i]=GetClusterMI(index);
2001   }
2002
2003   Int_t lab=123456789;
2004   for (i=0; i<noc; i++) {
2005     AliTPCclusterMI *c=clusters[i];
2006     if (!clusters[i]) continue;
2007     lab=TMath::Abs(c->GetLabel(0));
2008     Int_t j;
2009     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
2010     lb[j]=lab;
2011     (mx[j])++;
2012   }
2013
2014   Int_t max=0;
2015   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
2016
2017   for (i=0; i<noc; i++) {
2018     AliTPCclusterMI *c=clusters[i];
2019     if (!clusters[i]) continue;
2020     if (TMath::Abs(c->GetLabel(1)) == lab ||
2021         TMath::Abs(c->GetLabel(2)) == lab ) max++;
2022   }
2023
2024   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
2025
2026   else {
2027      Int_t tail=Int_t(0.10*noc);
2028      max=0;
2029      for (i=1; i<=tail; i++) {
2030        AliTPCclusterMI *c=clusters[noc-i];
2031        if (!clusters[i]) continue;
2032        if (lab == TMath::Abs(c->GetLabel(0)) ||
2033            lab == TMath::Abs(c->GetLabel(1)) ||
2034            lab == TMath::Abs(c->GetLabel(2))) max++;
2035      }
2036      if (max < Int_t(0.5*tail)) lab=-lab;
2037   }
2038
2039   t->SetLabel(lab);
2040
2041   delete[] lb;
2042   delete[] mx;
2043   delete[] clusters;
2044 }
2045
2046 //_________________________________________________________________________
2047 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
2048   //-----------------------------------------------------------------------
2049   // Setup inner sector
2050   //-----------------------------------------------------------------------
2051   if (f==0) {
2052      fAlpha=par->GetInnerAngle();
2053      fAlphaShift=par->GetInnerAngleShift();
2054      fPadPitchWidth=par->GetInnerPadPitchWidth();
2055      fPadPitchLength=par->GetInnerPadPitchLength();
2056      fN=par->GetNRowLow();
2057      fRow=new AliTPCRow[fN];
2058      for (Int_t i=0; i<fN; i++) {
2059        fRow[i].SetX(par->GetPadRowRadiiLow(i));
2060        fRow[i].fDeadZone =1.5;  //1.5 cm of dead zone
2061      }
2062   } else {
2063      fAlpha=par->GetOuterAngle();
2064      fAlphaShift=par->GetOuterAngleShift();
2065      fPadPitchWidth  = par->GetOuterPadPitchWidth();
2066      fPadPitchLength = par->GetOuter1PadPitchLength();
2067      f1PadPitchLength = par->GetOuter1PadPitchLength();
2068      f2PadPitchLength = par->GetOuter2PadPitchLength();
2069
2070      fN=par->GetNRowUp();
2071      fRow=new AliTPCRow[fN];
2072      for (Int_t i=0; i<fN; i++) {
2073        fRow[i].SetX(par->GetPadRowRadiiUp(i));
2074        fRow[i].fDeadZone =1.5;  // 1.5 cm of dead zone
2075      }
2076   }
2077 }
2078
2079
2080 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
2081   //
2082   if (fClusterTracks) delete [] fClusterTracks;
2083   fClusterTracks = 0;
2084 }
2085
2086 void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
2087   //create cluster tracks
2088   if (fN>0)
2089     fClusterTracks = new AliTPCclusterTracks[fN];
2090 }
2091
2092 void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
2093   if (fClusterTracks) delete[] fClusterTracks;
2094   fClusterTracks =0;
2095 }
2096
2097
2098
2099 void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
2100   //
2101   //
2102   // update information of the cluster tracks - if track is nearer then other tracks to the
2103   // given track
2104   const AliTPCclusterMI * cl = (*this)[clindex];
2105   AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
2106   // find the distance of the cluster to the track
2107   Float_t dy2 = (cl->GetY()- seed->GetY());
2108   dy2*=dy2;
2109   Float_t dz2 = (cl->GetZ()- seed->GetZ());
2110   dz2*=dz2;
2111   //
2112   Float_t distance = TMath::Sqrt(dy2+dz2);
2113   if (distance> 3)
2114     return;  // MI - to be changed - AliTPCtrackerParam
2115
2116   if ( distance < cltracks->fDistance[0]){
2117     cltracks->fDistance[2] =cltracks->fDistance[1];
2118     cltracks->fDistance[1] =cltracks->fDistance[0];
2119     cltracks->fDistance[0] =distance;
2120     cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
2121     cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
2122     cltracks->fTrackIndex[0] =trindex;
2123   }
2124   else
2125     if ( distance < cltracks->fDistance[1]){
2126       cltracks->fDistance[2] =cltracks->fDistance[1];
2127       cltracks->fDistance[1] =distance;
2128       cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
2129       cltracks->fTrackIndex[1] =trindex;
2130     } else
2131       if (distance < cltracks->fDistance[2]){
2132         cltracks->fDistance[2] =distance;
2133         cltracks->fTrackIndex[2] =trindex;
2134       }
2135 }
2136
2137
2138 //_________________________________________________________________________
2139 void
2140 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
2141   //-----------------------------------------------------------------------
2142   // Insert a cluster into this pad row in accordence with its y-coordinate
2143   //-----------------------------------------------------------------------
2144   if (fN==kMaxClusterPerRow) {
2145     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
2146   }
2147   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
2148   Int_t i=Find(c->GetZ());
2149   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
2150   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
2151   fIndex[i]=index; fClusters[i]=c; fN++;
2152 }
2153
2154 //___________________________________________________________________
2155 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
2156   //-----------------------------------------------------------------------
2157   // Return the index of the nearest cluster
2158   //-----------------------------------------------------------------------
2159   if (fN==0) return 0;
2160   if (z <= fClusters[0]->GetZ()) return 0;
2161   if (z > fClusters[fN-1]->GetZ()) return fN;
2162   Int_t b=0, e=fN-1, m=(b+e)/2;
2163   for (; b<e; m=(b+e)/2) {
2164     if (z > fClusters[m]->GetZ()) b=m+1;
2165     else e=m;
2166   }
2167   return m;
2168 }
2169
2170
2171
2172 //___________________________________________________________________
2173 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
2174   //-----------------------------------------------------------------------
2175   // Return the index of the nearest cluster in z y
2176   //-----------------------------------------------------------------------
2177   Float_t maxdistance = roady*roady + roadz*roadz;
2178
2179   AliTPCclusterMI *cl =0;
2180   for (Int_t i=Find(z-roadz); i<fN; i++) {
2181       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
2182       if (c->GetZ() > z+roadz) break;
2183       if ( (c->GetY()-y) >  roady ) continue;
2184       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
2185       if (maxdistance>distance) {
2186         maxdistance = distance;
2187         cl=c;
2188       }
2189   }
2190   return cl;
2191 }
2192
2193
2194
2195
2196
2197 AliTPCseed::AliTPCseed():AliTPCtrack(){
2198   //
2199   fRow=0;
2200   fRemoval =0;
2201   memset(fClusterIndex,0,sizeof(Int_t)*200);
2202   fPoints = 0;
2203   fEPoints = 0;
2204   fNFoundable =0;
2205   fNShared  =0;
2206   fTrackPoints =0;
2207   fRemoval = 0;
2208 }
2209
2210 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
2211   fPoints = 0;
2212   fEPoints = 0;
2213   fNShared  =0;
2214   fTrackPoints =0;
2215   fRemoval =0;
2216 }
2217
2218 AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
2219   fRow=0;
2220   memset(fClusterIndex,0,sizeof(Int_t)*200);
2221   fPoints = 0;
2222   fEPoints = 0;
2223   fNFoundable =0;
2224   fNShared  =0;
2225   fTrackPoints =0;
2226   fRemoval =0;
2227 }
2228
2229 AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
2230                                         Double_t xr, Double_t alpha):
2231   AliTPCtrack(index, xx, cc, xr, alpha) {
2232   //
2233   //
2234   fRow =0;
2235   memset(fClusterIndex,0,sizeof(Int_t)*200);
2236   fPoints = 0;
2237   fEPoints = 0;
2238   fNFoundable =0;
2239   fNShared  = 0;
2240   fTrackPoints =0;
2241   fRemoval =0;
2242 }
2243
2244 AliTPCseed::~AliTPCseed(){
2245   if (fPoints) delete fPoints;
2246   fPoints =0;
2247   fEPoints = 0;
2248   if (fTrackPoints){
2249     for (Int_t i=0;i<8;i++){
2250       delete [] fTrackPoints[i];
2251     }
2252     delete fTrackPoints;
2253     fTrackPoints =0;
2254   }
2255
2256 }
2257
2258 AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
2259 {
2260   //
2261   //
2262   if (!fTrackPoints) {
2263     fTrackPoints = new AliTPCTrackPoint*[8];
2264     for ( Int_t i=0;i<8;i++)
2265       fTrackPoints[i]=0;
2266   }
2267   Int_t index1 = i/20;
2268   if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
2269   return &(fTrackPoints[index1][i%20]);
2270 }
2271
2272 void AliTPCseed::RebuildSeed()
2273 {
2274   //
2275   // rebuild seed to be ready for storing
2276   fPoints = new TClonesArray("AliTPCTrackPoint",160);
2277   fPoints->ExpandCreateFast(160);
2278   fEPoints = new TClonesArray("AliTPCExactPoint",1);
2279   for (Int_t i=0;i<160;i++){
2280     AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
2281     *trpoint = *(GetTrackPoint(i));
2282   }
2283
2284 }
2285
2286 //_____________________________________________________________________________
2287 void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
2288   //-----------------------------------------------------------------
2289   // This funtion calculates dE/dX within the "low" and "up" cuts.
2290   //-----------------------------------------------------------------
2291
2292   Float_t amp[200];
2293   Float_t angular[200];
2294   Float_t weight[200];
2295   Int_t index[200];
2296   //Int_t nc = 0;
2297   //  TClonesArray & arr = *fPoints;
2298   Float_t meanlog = 100.;
2299
2300   Float_t mean[4]  = {0,0,0,0};
2301   Float_t sigma[4] = {1000,1000,1000,1000};
2302   Int_t nc[4]      = {0,0,0,0};
2303   Float_t norm[4]    = {1000,1000,1000,1000};
2304   //
2305   //
2306   fNShared =0;
2307
2308   for (Int_t of =0; of<4; of++){
2309     for (Int_t i=of;i<160;i+=4)
2310       {
2311         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
2312         AliTPCTrackPoint * point = GetTrackPoint(i);
2313         if (point==0) continue;
2314         if (point->fIsShared){
2315           fNShared++;
2316           continue;
2317         }
2318         if (point->GetCPoint().GetMax()<5) continue;
2319         Float_t angley = point->GetTPoint().GetAngleY();
2320         Float_t anglez = point->GetTPoint().GetAngleZ();
2321         Int_t   type   = point->GetCPoint().GetType();
2322         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
2323         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
2324         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
2325
2326         Float_t ampc   = 0;     // normalization to the number of electrons
2327         if (i>64){
2328           ampc = 1.*point->GetCPoint().GetMax();
2329           //ampc = 1.*point->GetCPoint().GetQ();
2330           //      AliTPCClusterPoint & p = point->GetCPoint();
2331           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
2332           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2333           //Float_t dz =
2334           //  TMath::Abs( Int_t(iz) - iz + 0.5);
2335           //ampc *= 1.15*(1-0.3*dy);
2336           //ampc *= 1.15*(1-0.3*dz);
2337           //      Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
2338           //ampc               *=zfactor;
2339         }
2340         else{
2341           ampc = 1.0*point->GetCPoint().GetMax();
2342           //ampc = 1.0*point->GetCPoint().GetQ();
2343           //AliTPCClusterPoint & p = point->GetCPoint();
2344           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
2345           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2346           //Float_t dz =
2347           //  TMath::Abs( Int_t(iz) - iz + 0.5);
2348
2349           //ampc *= 1.15*(1-0.3*dy);
2350           //ampc *= 1.15*(1-0.3*dz);
2351           //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
2352           //ampc               *=zfactor;
2353
2354         }
2355         ampc *= 2.0;     // put mean value to channel 50
2356         //ampc *= 0.58;     // put mean value to channel 50
2357         Float_t w      =  1.;
2358         //      if (type>0)  w =  1./(type/2.-0.5);
2359         //      Float_t z = TMath::Abs(point->GetCPoint().GetZ());
2360         if (i<64) {
2361           ampc /= 0.6;
2362           //ampc /= (1+0.0008*z);
2363         } else
2364           if (i>128){
2365             ampc /=1.5;
2366             //ampc /= (1+0.0008*z);
2367           }else{
2368             //ampc /= (1+0.0008*z);
2369           }
2370
2371         if (type<0) {  //amp at the border - lower weight
2372           // w*= 2.;
2373
2374           continue;
2375         }
2376         if (rsigma>1.5) ampc/=1.3;  // if big backround
2377         amp[nc[of]]        = ampc;
2378         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
2379         weight[nc[of]]     = w;
2380         nc[of]++;
2381       }
2382
2383     TMath::Sort(nc[of],amp,index,kFALSE);
2384     Float_t sumamp=0;
2385     Float_t sumamp2=0;
2386     Float_t sumw=0;
2387     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
2388     meanlog = 200;
2389     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
2390       Float_t ampl      = amp[index[i]]/angular[index[i]];
2391       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
2392       //
2393       sumw    += weight[index[i]];
2394       sumamp  += weight[index[i]]*ampl;
2395       sumamp2 += weight[index[i]]*ampl*ampl;
2396       norm[of]    += angular[index[i]]*weight[index[i]];
2397     }
2398     if (sumw<1){
2399       SetdEdx(0);
2400     }
2401     else {
2402       norm[of] /= sumw;
2403       mean[of]  = sumamp/sumw;
2404       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
2405       if (sigma[of]>0.1)
2406         sigma[of] = TMath::Sqrt(sigma[of]);
2407       else
2408         sigma[of] = 1000;
2409
2410     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
2411     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
2412     //mean *=(1-0.1*(norm-1.));
2413     }
2414   }
2415
2416   Float_t dedx =0;
2417   fSdEdx =0;
2418   fMAngular =0;
2419   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
2420   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
2421
2422
2423   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
2424   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
2425
2426   Int_t norm2 = 0;
2427   Int_t norm3 = 0;
2428   for (Int_t i =0;i<4;i++){
2429     if (nc[i]>2&&nc[i]<1000){
2430       dedx      += mean[i] *nc[i];
2431       fSdEdx    += sigma[i]*(nc[i]-2);
2432       fMAngular += norm[i] *nc[i];
2433       norm2     += nc[i];
2434       norm3     += nc[i]-2;
2435     }
2436     fDEDX[i]  = mean[i];
2437     fSDEDX[i] = sigma[i];
2438     fNCDEDX[i]= nc[i];
2439   }
2440
2441   if (norm3>0){
2442     dedx   /=norm2;
2443     fSdEdx /=norm3;
2444     fMAngular/=norm2;
2445   }
2446   else{
2447     SetdEdx(0);
2448     return;
2449   }
2450   //  Float_t dedx1 =dedx;
2451
2452   dedx =0;
2453   for (Int_t i =0;i<4;i++){
2454     if (nc[i]>2&&nc[i]<1000){
2455       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
2456       dedx      += mean[i] *nc[i];
2457     }
2458     fDEDX[i]  = mean[i];
2459   }
2460   dedx /= norm2;
2461
2462
2463
2464   SetdEdx(dedx);
2465
2466   //mi deDX
2467
2468
2469
2470   //Very rough PID
2471   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
2472
2473   if (p<0.6) {
2474     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2475     if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
2476     SetMass(0.93827); return;
2477   }
2478
2479   if (p<1.2) {
2480     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2481     SetMass(0.93827); return;
2482   }
2483
2484   SetMass(0.13957); return;
2485
2486 }
2487
2488