TPC/AliTPCtrackerMI.cxx

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