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