TPC/AliTPCtrackerMI.cxx

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