TPC/AliTPCtrackerMI.cxx

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