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