TPC/AliTPCtrackerMI.cxx

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