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   Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
 538   for (Int_t i=0; i<j; i++) {
 539     //  AliSegmentID *s=fClustersArray.LoadEntry(i);
 540     AliSegmentID *s= const_cast<AliSegmentID*>(fClustersArray.At(i));
 541     if (!s) continue;
 542     Int_t sec,row;
 543     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 544     par->AdjustSectorRow(s->GetID(),sec,row);
 545     if (sec<fkNIS*2) continue;
 546     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 547     Int_t ncl=clrow->GetArray()->GetEntriesFast();
 548     while (ncl--) {
 549       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
 550       index=(((sec<<8)+row)<<16)+ncl;
 551       fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
 552     }
 553   }
 554   fN=fkNOS;
 555   fSectors=fOuterSec;
 556 }
 557
 558
 559 //_____________________________________________________________________________
 560 void AliTPCtrackerMI::LoadInnerSectors() {
 561   //-----------------------------------------------------------------
 562   // This function fills inner TPC sectors with clusters.
 563   //-----------------------------------------------------------------
 564   UInt_t index;
 565   //Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
 566   Int_t j = ((AliTPCParam*)fParam)->GetNRowsTotal();
 567   for (Int_t i=0; i<j; i++) {
 568     //   AliSegmentID *s=fClustersArray.LoadEntry(i);
 569     AliSegmentID *s=const_cast<AliSegmentID*>(fClustersArray.At(i));
 570     if (!s) continue;
 571     Int_t sec,row;
 572     AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
 573     par->AdjustSectorRow(s->GetID(),sec,row);
 574     if (sec>=fkNIS*2) continue;
 575     AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
 576     Int_t ncl=clrow->GetArray()->GetEntriesFast();
 577     while (ncl--) {
 578       AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
 579       index=(((sec<<8)+row)<<16)+ncl;
 580       fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
 581     }
 582   }
 583
 584   fN=fkNIS;
 585   fSectors=fInnerSec;
 586 }
 587
 588 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
 589   //-----------------------------------------------------------------
 590   // This function tries to find a track prolongation to next pad row
 591   //-----------------------------------------------------------------
 592   //  Double_t xt=t.GetX();
 593   //  Int_t row = fSectors->GetRowNumber(xt)-1;
 594   //  if (row < nr) return 1; // don't prolongate if not information until now -
 595   //
 596   Double_t  x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
 597   //  if (t.GetRadius()>x+10 ) return 0;
 598
 599   if (!t.PropagateTo(x)) {
 600     t.fStopped = kTRUE;
 601     return 0;
 602   }
 603   // update current
 604   t.fCurrentSigmaY = GetSigmaY(&t);
 605   t.fCurrentSigmaZ = GetSigmaZ(&t);
 606   //
 607   AliTPCclusterMI *cl=0;
 608   UInt_t index=0;
 609   const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
 610   Double_t sy2=ErrY2(&t)*2;
 611   Double_t sz2=ErrZ2(&t)*2;
 612
 613
 614   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
 615   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
 616   Double_t  y=t.GetY(), z=t.GetZ();
 617
 618   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
 619     t.fInDead = kTRUE;
 620     return 0;
 621   }
 622   else
 623     {
 624       if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
 625       else
 626         return 0;
 627     }
 628   //calculate
 629   Float_t maxdistance = roady*roady + roadz*roadz;
 630   if (krow) {
 631     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
 632       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
 633       if (c->GetZ() > z+roadz) break;
 634       if ( (c->GetY()-y) >  roady ) continue;
 635       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
 636       if (maxdistance>distance) {
 637         maxdistance = distance;
 638         cl=c;
 639         index=krow.GetIndex(i);
 640       }
 641     }
 642   }
 643   if (cl) {
 644     //Double_t sy2=
 645     ErrY2(&t,cl);
 646     //Double_t sz2=
 647     ErrZ2(&t,cl);
 648     Double_t chi2= t.GetPredictedChi2(cl);
 649     UpdateTrack(&t,cl,chi2,index);
 650   } else {
 651     if (y > ymax) {
 652       t.fRelativeSector= (t.fRelativeSector+1) % fN;
 653       if (!t.Rotate(fSectors->GetAlpha()))
 654         return 0;
 655     } else if (y <-ymax) {
 656       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
 657       if (!t.Rotate(-fSectors->GetAlpha()))
 658         return 0;
 659     }
 660   }
 661   return 1;
 662 }
 663
 664
 665 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex,  Int_t nr) {
 666   //-----------------------------------------------------------------
 667   // This function tries to find a track prolongation to next pad row
 668   //-----------------------------------------------------------------
 669   t.fCurrentCluster  = 0;
 670   t.fCurrentClusterIndex1 = 0;
 671   t.fCurrentClusterIndex2 = 0;
 672
 673   Double_t xt=t.GetX();
 674   Int_t row = fSectors->GetRowNumber(xt)-1;
 675   if (row < nr) return 1; // don't prolongate if not information until now -
 676   Double_t x=fSectors->GetX(nr);
 677   //  if (t.fStopped) return 0;
 678   //  if (t.GetRadius()>x+10 ) return 0;
 679   if (!t.PropagateTo(x)){
 680     t.fStopped =kTRUE;
 681     return 0;
 682   }
 683   // update current
 684   t.fCurrentSigmaY = GetSigmaY(&t);
 685   t.fCurrentSigmaZ = GetSigmaZ(&t);
 686
 687   AliTPCclusterMI *cl=0;
 688   UInt_t index=0;
 689   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
 690   //
 691   Double_t  y=t.GetY(), z=t.GetZ();
 692   Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
 693   Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
 694   //
 695
 696   Float_t maxdistance = 1000000;
 697   if (krow) {
 698     for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
 699       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
 700       if (c->GetZ() > z+roadz) break;
 701       if (TMath::Abs(c->GetY()-y)>roady) continue;
 702
 703       //krow.UpdateClusterTrack(i,trindex,&t);
 704
 705       Float_t dy2 = (c->GetY()- t.GetY());
 706       dy2*=dy2;
 707       Float_t dz2 = (c->GetZ()- t.GetZ());
 708       dz2*=dz2;
 709       //
 710       Float_t distance = dy2+dz2;
 711       //
 712       if (distance > maxdistance) continue;
 713       maxdistance = distance;
 714       cl=c;
 715       index=i;
 716     }
 717   }
 718   t.fCurrentCluster  = cl;
 719   t.fCurrentClusterIndex1 = krow.GetIndex(index);
 720   t.fCurrentClusterIndex2 = index;
 721   return 1;
 722 }
 723
 724
 725 Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
 726   //-----------------------------------------------------------------
 727   // This function tries to find a track prolongation to next pad row
 728   //-----------------------------------------------------------------
 729   AliTPCseed & t  = *((AliTPCseed*)(fSeeds->At(trindex)));
 730   AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
 731   //  Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
 732   Double_t y=t.GetY();
 733   Double_t ymax=fSectors->GetMaxY(nr);
 734
 735   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
 736     t.fInDead = kTRUE;
 737     return 0;
 738   }
 739   else
 740     {
 741       if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
 742       else
 743         return 0;
 744     }
 745
 746   if (t.fCurrentCluster) {
 747     //    Float_t l=fSectors->GetPadPitchLength();
 748     //    AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
 749
 750     Double_t sy2=ErrY2(&t,t.fCurrentCluster);
 751     Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
 752
 753
 754     Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
 755     Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
 756
 757     Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
 758     Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
 759
 760     Double_t rdistance  = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
 761
 762
 763     //    printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
 764     if ( (rdistancey>1) || (rdistancez>1)) return 0;
 765     if (rdistance>4) return 0;
 766
 767     if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)
 768         return 0;  //suspisiouce - will be changed
 769
 770     if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)
 771         // strict cut on overlaped cluster
 772         return 0;  //suspisiouce - will be changed
 773
 774     if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 )
 775          && t.fCurrentCluster->GetType()<0)
 776       return 0;
 777
 778     //    t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
 779     UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
 780
 781   } else {
 782     if (y > ymax) {
 783       t.fRelativeSector= (t.fRelativeSector+1) % fN;
 784       if (!t.Rotate(fSectors->GetAlpha()))
 785         return 0;
 786     } else if (y <-ymax) {
 787       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
 788       if (!t.Rotate(-fSectors->GetAlpha()))
 789         return 0;
 790     }
 791   }
 792   return 1;
 793 }
 794
 795
 796
 797
 798 /*
 799 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
 800 {
 801   //-----------------------------------------------------------------
 802   // fast prolongation mathod -
 803   // don't update track only after step clusters
 804   //-----------------------------------------------------------------
 805   Double_t xt=t.GetX();
 806   //
 807   Double_t alpha=t.GetAlpha();
 808   alpha =- fSectors->GetAlphaShift();
 809   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 810   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 811   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
 812   Int_t row0 = fSectors->GetRowNumber(xt);
 813   Double_t x    = fSectors->GetX(row0);
 814   Double_t ymax = fSectors->GetMaxY(row0);
 815   //
 816   Double_t sy2=ErrY2(&t)*2;
 817   Double_t sz2=ErrZ2(&t)*2;
 818   Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
 819   Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
 820   Float_t maxdistance = roady*roady + roadz*roadz;
 821   t.fCurrentSigmaY = GetSigmaY(&t);
 822   t.fCurrentSigmaZ = GetSigmaZ(&t);
 823   //
 824   Int_t nclusters = 0;
 825   Double_t y;
 826   Double_t z;
 827   Double_t yy[200];   //track prolongation
 828   Double_t zz[200];
 829   Double_t cy[200];  // founded cluster position
 830   Double_t cz[200];
 831   Double_t sy[200];  // founded cluster error
 832   Double_t sz[200];
 833   Bool_t   hitted[200]; // indication of cluster presence
 834   //
 835
 836   //
 837   for (Int_t drow = step; drow>=0; drow--) {
 838     Int_t row = row0-drow;
 839     if (row<0) break;
 840     Double_t x    = fSectors->GetX(row);
 841     Double_t ymax = fSectors->GetMaxY(row);
 842     t.GetProlongation(x,y,z);
 843     yy[drow] =y;
 844     zz[drow] =z;
 845     const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
 846     if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
 847       t.fInDead = kTRUE;
 848       break;
 849     }
 850     else
 851       {
 852         t.fNFoundable++;
 853       }
 854
 855     //find nearest  cluster
 856     AliTPCclusterMI *cl= 0;
 857     if (krow) {
 858       for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
 859         AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
 860         if (c->GetZ() > z+roadz) break;
 861         if ( (c->GetY()-y) >  roady ) continue;
 862         Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
 863         if (maxdistance>distance) {
 864           maxdistance = distance;
 865           cl=c;
 866           //      index=krow.GetIndex(i);
 867         }
 868       }
 869     }  // end of seearch
 870     //update cluster information
 871     if (cl){
 872       cy[drow] = cl->GetY();
 873       cz[drow] = cl->GetZ();
 874       sy[drow] = ErrY2(&t,cl);
 875       sz[drow] = ErrZ2(&t,cl);
 876       hitted[drow] = kTRUE;
 877       nclusters++;
 878     }
 879     else
 880       hitted[drow] = kFALSE;
 881   }
 882   //if we have information - update track
 883   if (nclusters>0){
 884     Float_t sumyw0   = 0;
 885     Float_t sumzw0   = 0;
 886     Float_t sumyw   = 0;
 887     Float_t sumzw   = 0;
 888     Float_t sumrow  = 0;
 889     for (Int_t i=0;i<step;i++)
 890       if (hitted[i]){
 891         sumyw0+= 1/sy[i];
 892         sumzw0+= 1/sz[i];
 893         //
 894         sumyw+= (cy[i]-yy[i])/sy[i];
 895         sumzw+= (cz[i]-zz[i])/sz[i];
 896         sumrow+= i;
 897       }
 898     Float_t dy = sumyw/sumyw0;
 899     Float_t dz = sumzw/sumzw0;
 900     Float_t mrow = sumrow/nclusters+row0;
 901     Float_t x = fSectors->GetX(mrow);
 902     t.PropagateTo(x);
 903     AliTPCclusterMI cvirtual;
 904     cvirtual.SetZ(dz+t.GetZ());
 905     cvirtual.SetY(dy+t.GetY());
 906     t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
 907     t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
 908     Float_t chi2 = t.GetPredictedChi2(&cvirtual);
 909     t.Update(&cvirtual,chi2,0);
 910     Int_t ncl = t.GetNumberOfClusters();
 911     ncl = ncl-1+nclusters;
 912     t.SetN(ncl);
 913   }
 914   return  1;
 915 }
 916 */
 917
 918 //_____________________________________________________________________________
 919 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
 920   //-----------------------------------------------------------------
 921   // This function tries to find a track prolongation.
 922   //-----------------------------------------------------------------
 923   Double_t xt=t.GetX();
 924   //
 925   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
 926   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 927   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 928   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
 929
 930   for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
 931
 932     if (FollowToNext(t,nr)==0) {
 933     }
 934   }
 935   return 1;
 936 }
 937
 938
 939 //_____________________________________________________________________________
 940 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
 941   //-----------------------------------------------------------------
 942   // This function tries to find a track prolongation.
 943   //-----------------------------------------------------------------
 944   Double_t xt=t.GetX();
 945   //
 946   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
 947   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
 948   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
 949   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
 950
 951   for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
 952     if (t.GetSnp()<0.8)
 953       FollowToNext(t,nr);
 954   }
 955   return 1;
 956 }
 957
 958
 959
 960
 961
 962 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
 963 {
 964   //
 965   //
 966   sum1=0;
 967   sum2=0;
 968   Int_t sum=0;
 969   if (s1->fSector!=s2->fSector) return 0;
 970   //
 971   Float_t dz2 =(s1->GetZ() - s2->GetZ());
 972   dz2*=dz2;
 973   Float_t dy2 =(s1->GetY() - s2->GetY());
 974   dy2*=dy2;
 975   Float_t distance = TMath::Sqrt(dz2+dy2);
 976   if (distance>5.) return 0; // if there are far away  - not overlap - to reduce combinatorics
 977
 978   for (Int_t i=0;i<160;i++){
 979     if (s1->fClusterIndex[i]>0) sum1++;
 980     if (s2->fClusterIndex[i]>0) sum2++;
 981     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
 982       sum++;
 983     }
 984   }
 985
 986   Float_t summin = TMath::Min(sum1+1,sum2+1);
 987   Float_t ratio = (sum+1)/Float_t(summin);
 988   return ratio;
 989 }
 990
 991 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
 992 {
 993   //
 994   //
 995   if (s1->fSector!=s2->fSector) return;
 996   //
 997   Float_t dz2 =(s1->GetZ() - s2->GetZ());
 998   dz2*=dz2;
 999   Float_t dy2 =(s1->GetY() - s2->GetY());
1000   dy2*=dy2;
1001   Float_t distance = TMath::Sqrt(dz2+dy2);
1002   if (distance>15.) return ; // if there are far away  - not overlap - to reduce combinatorics
1003   //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
1004   //  TClonesArray &pointarray1 = *(s1->fPoints);
1005   //TClonesArray &pointarray2 = *(s2->fPoints);
1006   //
1007   for (Int_t i=0;i<160;i++){
1008     if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
1009       //  AliTPCTrackPoint *p1  = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
1010       //AliTPCTrackPoint *p2  = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i));
1011       AliTPCTrackPoint *p1  = s1->GetTrackPoint(i);
1012       AliTPCTrackPoint *p2  = s2->GetTrackPoint(i);;
1013       p1->fIsShared = kTRUE;
1014       p2->fIsShared = kTRUE;
1015     }
1016   }
1017 }
1018
1019
1020
1021
1022 void  AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared){
1023
1024
1025
1026   //
1027   // remove overlap - used removal factor - removal index stored in the track
1028   arr->Sort();  // sorting according z
1029   arr->Expand(arr->GetEntries());
1030   Int_t nseed=arr->GetEntriesFast();
1031   //  printf("seeds \t%p \t%d\n",arr, nseed);
1032   //  arr->Expand(arr->GetEntries());  //remove 0 pointers
1033   nseed = arr->GetEntriesFast();
1034   Int_t removed = 0;
1035   for (Int_t i=0; i<nseed; i++) {
1036     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1037     if (!pt) {
1038       continue;
1039     }
1040     if (!(pt->IsActive())) continue;
1041     for (Int_t j=i+1; j<nseed; j++){
1042       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1043       if ((pt2) && pt2->IsActive())
1044         if (pt->fSector == pt2->fSector)
1045           if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
1046             Int_t sum1,sum2;
1047             Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
1048             //if (sum1==0) {
1049             //  pt->Desactivate(removalindex); // arr->RemoveAt(i);
1050             //  break;
1051             //}
1052             if (ratio>factor){
1053               //          if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex);  // arr->RemoveAt(j);
1054               Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
1055               Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
1056               removed++;
1057               if (TMath::Abs(ratio3)>0.025){  // if much more points
1058                 if (sum1>sum2) pt2->Desactivate(removalindex);
1059                 else {
1060                   pt->Desactivate(removalindex); // arr->RemoveAt(i);
1061                   break;
1062                 }
1063               }
1064               else{  //decide on mean chi2
1065                 if (ratio2<1)
1066                   pt2->Desactivate(removalindex);
1067                 else {
1068                   pt->Desactivate(removalindex); // arr->RemoveAt(i);
1069                   break;
1070                 }
1071               }
1072
1073             }  // if suspicious ratio
1074           }
1075           else
1076             break;
1077     }
1078   }
1079   //  printf("removed\t%d\n",removed);
1080   Int_t good =0;
1081   for (Int_t i=0; i<nseed; i++) {
1082     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1083     if (!pt) break;
1084     if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
1085       //desactivate tracks with small number of points
1086       //      printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
1087       pt->Desactivate(10);  //desactivate  - small muber of points
1088     }
1089     if (!(pt->IsActive())) continue;
1090     good++;
1091   }
1092
1093
1094   if (shared)
1095     for (Int_t i=0; i<nseed; i++) {
1096       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1097       if (!pt) continue;
1098       if (!(pt->IsActive())) continue;
1099       for (Int_t j=i+1; j<nseed; j++){
1100         AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1101         if ((pt2) && pt2->IsActive()) {
1102           if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
1103           SignShared(pt,pt2);
1104         }
1105       }
1106     }
1107   fNtracks = good;
1108   printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
1109
1110
1111 }
1112 void AliTPCtrackerMI::MakeSeedsAll()
1113 {
1114   if (fSeeds == 0) fSeeds = new TObjArray;
1115   TObjArray * arr;
1116   for (Int_t sec=0;sec<fkNOS;sec+=3){
1117      arr = MakeSeedsSectors(sec,sec+3);
1118      Int_t nseed = arr->GetEntriesFast();
1119      for (Int_t i=0;i<nseed;i++)
1120        fSeeds->AddLast(arr->RemoveAt(i));
1121   }
1122   //  fSeeds = MakeSeedsSectors(0,fkNOS);
1123 }
1124
1125 TObjArray *  AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
1126 {
1127   //
1128   // loop over all  sectors and make seed
1129   //find track seeds
1130   TStopwatch timer;
1131   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1132   Int_t nrows=nlow+nup;
1133   Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
1134   //  if (fSeeds==0) fSeeds = new TObjArray;
1135   TObjArray * arr = new TObjArray;
1136
1137   for (Int_t sec=sec1; sec<sec2;sec++){
1138     MakeSeeds(arr, sec, nup-1, nup-1-gap);
1139     MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
1140   }
1141   gap = Int_t(0.3* nrows);
1142   for (Int_t sec=sec1; sec<sec2;sec++){
1143     //find secondaries
1144     MakeSeeds2(arr, sec, nup-1, nup-1-gap);
1145     MakeSeeds2(arr, sec, nup-1-shift, nup-1-shift-gap);
1146     MakeSeeds2(arr, sec, nup-1-2*shift, nup-1-2*shift-gap);
1147     //MakeSeeds2(arr, sec, nup-1-3*shift, nup-1-3*shift-gap);
1148     MakeSeeds2(arr, sec, 30, 0);
1149   }
1150
1151   Int_t nseed=arr->GetEntriesFast();
1152   gap=Int_t(0.3*nrows);
1153   // continue seeds
1154   Int_t i;
1155   for (i=0; i<nseed; i++) {
1156     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
1157     if (!pt) continue;
1158     if (FollowProlongation(t,nup-gap)) {
1159       pt->fIsSeeding =kFALSE;
1160       continue;
1161     }
1162     delete arr->RemoveAt(i);
1163   }
1164
1165   //
1166   //remove seeds which overlaps
1167   RemoveOverlap(arr,0.6,1);
1168   //delete seeds - which were sign
1169   nseed=arr->GetEntriesFast();
1170   for (i=0; i<nseed; i++) {
1171     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1172     //, &t=*pt;
1173     if (!pt) continue;
1174     if ((pt->IsActive())  &&  pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
1175       //pt->Reset();
1176       //FollowBackProlongation(*pt,nup-1);
1177       //if ( pt->GetNumberOfClusters() < pt->fNFoundable*0.5 || pt->GetNumberOfClusters()<10 )
1178       //delete arr->RemoveAt(i);
1179       //else
1180       //        pt->Reset();
1181       continue;
1182     }
1183     delete arr->RemoveAt(i);
1184   }
1185   //RemoveOverlap(arr,0.6,1);
1186   return arr;
1187 }
1188
1189
1190
1191 //_____________________________________________________________________________
1192 void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1193   //-----------------------------------------------------------------
1194   // This function creates track seeds.
1195   //-----------------------------------------------------------------
1196   //  if (fSeeds==0) fSeeds=new TObjArray(15000);
1197
1198   Double_t x[5], c[15];
1199
1200   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1201   Double_t cs=cos(alpha), sn=sin(alpha);
1202
1203   Double_t x1 =fOuterSec->GetX(i1);
1204   Double_t xx2=fOuterSec->GetX(i2);
1205   //
1206   //  for (Int_t ns=0; ns<fkNOS; ns++)
1207   Int_t ns =sec;
1208     {
1209     Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1210     Int_t nm=fOuterSec[ns][i2];
1211     Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
1212     const AliTPCRow& kr1=fOuterSec[ns][i1];
1213     AliTPCRow&  kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
1214     AliTPCRow&  kr22 = fOuterSec[(ns)%fkNOS][i2];
1215     AliTPCRow&  kr23 = fOuterSec[(ns+1)%fkNOS][i2];
1216
1217     for (Int_t is=0; is < kr1; is++) {
1218       Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
1219       Double_t x3=GetX(), y3=GetY(), z3=GetZ();
1220
1221       Float_t anglez = (z1-z3)/(x1-x3);
1222       Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
1223
1224       for (Int_t js=0; js < nl+nm+nu; js++) {
1225         const AliTPCclusterMI *kcl;
1226         Double_t x2,   y2,   z2;
1227         if (js<nl) {
1228           if (js==0) {
1229             js = kr21.Find(extraz-15.);
1230             if (js>=nl) continue;
1231           }
1232           kcl=kr21[js];
1233           z2=kcl->GetZ();
1234           if ((extraz-z2)>10) continue;
1235           if ((extraz-z2)<-10) {
1236             js = nl-1;
1237             continue;
1238           }
1239           y2=kcl->GetY();
1240           x2= xx2*cs+y2*sn;
1241           y2=-xx2*sn+y2*cs;
1242         } else
1243           if (js<nl+nm) {
1244             if (js==nl) {
1245               js = nl+kr22.Find(extraz-15.);
1246               if (js>=nl+nm) continue;
1247             }
1248             kcl=kr22[js-nl];
1249             z2=kcl->GetZ();
1250             if ((extraz-z2)>10) continue;
1251             if ((extraz-z2)<-10) {
1252               js = nl+nm-1;
1253               continue;
1254             }
1255             x2=xx2; y2=kcl->GetY();
1256           } else {
1257             //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];
1258             if (js==nl+nm) {
1259               js = nl+nm+kr23.Find(extraz-15.);
1260               if (js>=nl+nm+nu) break;
1261             }
1262             kcl=kr23[js-nl-nm];
1263             z2=kcl->GetZ();
1264             if ((extraz-z2)>10) continue;
1265             if ((extraz-z2)<-10) {
1266               break;
1267             }
1268             y2=kcl->GetY();
1269             x2=xx2*cs-y2*sn;
1270             y2=xx2*sn+y2*cs;
1271           }
1272
1273         Double_t zz=z1 - anglez*(x1-x2);
1274         if (TMath::Abs(zz-z2)>10.) continue;
1275
1276         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
1277         if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
1278
1279         x[0]=y1;
1280         x[1]=z1;
1281         x[4]=f1(x1,y1,x2,y2,x3,y3);
1282         if (TMath::Abs(x[4]) >= 0.0066) continue;
1283         x[2]=f2(x1,y1,x2,y2,x3,y3);
1284         //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1285         x[3]=f3(x1,y1,x2,y2,z1,z2);
1286         if (TMath::Abs(x[3]) > 1.2) continue;
1287         Double_t a=asin(x[2]);
1288         Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
1289         if (TMath::Abs(zv-z3)>10.) continue;
1290
1291         Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
1292         Double_t sy2=kcl->GetSigmaY2()*2,     sz2=kcl->GetSigmaZ2()*4;
1293         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1294         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1295         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1296
1297         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1298         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1299         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1300         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1301         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1302         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1303         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1304         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1305         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1306         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1307
1308         c[0]=sy1;
1309         c[1]=0.;       c[2]=sz1;
1310         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1311         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1312                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1313         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1314         c[13]=f30*sy1*f40+f32*sy2*f42;
1315         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1316
1317         UInt_t index=kr1.GetIndex(is);
1318         AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
1319         track->fIsSeeding = kTRUE;
1320         Int_t rc=FollowProlongation(*track, i2);
1321         //FollowProlongationFast(*track, 5);
1322         //FollowProlongationFast(*track, 5);
1323         //FollowProlongationFast(*track, 5);
1324         //FollowProlongationFast(*track, 5);
1325         //Int_t rc = 1;
1326
1327         track->fLastPoint = i1;  // first cluster in track position
1328         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1329         else arr->AddLast(track);
1330       }
1331     }
1332   }
1333 }
1334
1335
1336 //_____________________________________________________________________________
1337 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
1338   //-----------------------------------------------------------------
1339   // This function creates track seeds - without vertex constraint
1340   //-----------------------------------------------------------------
1341
1342   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
1343   //  Double_t cs=cos(alpha), sn=sin(alpha);
1344   Int_t row0 = (i1+i2)/2;
1345   Int_t drow = (i1-i2)/2;
1346   const AliTPCRow& kr0=fSectors[sec][row0];
1347   const AliTPCRow& krm=fSectors[sec][row0-1];
1348   const AliTPCRow& krp=fSectors[sec][row0+1];
1349   AliTPCRow * kr=0;
1350
1351   AliTPCpolyTrack polytrack;
1352   Int_t nclusters=fSectors[sec][row0];
1353
1354   for (Int_t is=0; is < nclusters; is++) {
1355     const AliTPCclusterMI * cl= kr0[is];
1356     Double_t x = kr0.GetX();
1357
1358     // Initialization of the polytrack
1359     polytrack.Reset();
1360
1361     Double_t y0= cl->GetY();
1362     Double_t z0= cl->GetZ();
1363     polytrack.AddPoint(x,y0,z0);
1364     Float_t roady = 5*TMath::Sqrt(cl->GetSigmaY2()+0.2);
1365     Float_t roadz = 5*TMath::Sqrt(cl->GetSigmaZ2()+0.2);
1366     //
1367     x = krm.GetX();
1368     cl = krm.FindNearest(y0,z0,roady,roadz);
1369     if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),roady,roadz);
1370     //
1371     x = krp.GetX();
1372     cl = krp.FindNearest(y0,z0,roady,roadz);
1373     if (cl) polytrack.AddPoint(x,cl->GetY(),cl->GetZ(),cl->GetSigmaY2()+0.05,cl->GetSigmaZ2()+0.05);
1374     //
1375     polytrack.UpdateParameters();
1376     // follow polytrack
1377     roadz = 0.6;
1378     roady = 0.6;
1379     //
1380     Double_t yn,zn;
1381     Int_t nfoundable = polytrack.GetN();
1382     Int_t nfound     = nfoundable;
1383     for (Int_t ddrow = 2; ddrow<drow;ddrow++){
1384       for (Int_t delta = -1;delta<=1;delta+=2){
1385         Int_t row = row0+ddrow*delta;
1386         kr = &(fSectors[sec][row]);
1387         Double_t xn = kr->GetX();
1388         Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone;
1389         polytrack.GetFitPoint(xn,yn,zn);
1390         if (TMath::Abs(yn)>ymax) continue;
1391         nfoundable++;
1392         AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
1393         if (cln) {
1394           polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),cln->GetSigmaY2()+0.05,cln->GetSigmaZ2()+0.05);
1395           nfound++;
1396         }
1397       }
1398       polytrack.UpdateParameters();
1399     }
1400     if ((nfound>0.5*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
1401       // add polytrack candidate
1402       Double_t x[5], c[15];
1403       Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
1404       polytrack.GetBoundaries(x3,x1);
1405       x2 = (x1+x3)/2.;
1406       polytrack.GetFitPoint(x1,y1,z1);
1407       polytrack.GetFitPoint(x2,y2,z2);
1408       polytrack.GetFitPoint(x3,y3,z3);
1409       //
1410       //is track pointing to the vertex ?
1411       Double_t x0,y0,z0;
1412       x0=0;
1413       polytrack.GetFitPoint(x0,y0,z0);
1414       if ( (TMath::Abs(z0-GetZ())<10) && (TMath::Abs(y0-GetY())<5)){ //if yes apply vertex constraint
1415         x3 = 0;
1416         y3 = GetY();
1417         z3 = GetZ();
1418       }
1419
1420       x[0]=y1;
1421       x[1]=z1;
1422       x[4]=f1(x1,y1,x2,y2,x3,y3);
1423       if (TMath::Abs(x[4]) >= 0.0066) continue;
1424       x[2]=f2(x1,y1,x2,y2,x3,y3);
1425       //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
1426       x[3]=f3(x1,y1,x2,y2,z1,z2);
1427       if (TMath::Abs(x[3]) > 1.2) continue;
1428       if (TMath::Abs(x[2]) > 0.99) continue;
1429       //      Double_t a=asin(x[2]);
1430
1431
1432       Double_t sy=1.5, sz=1.5;
1433       Double_t sy1=1.5, sz1=1.5;
1434       Double_t sy2=1.3, sz2=1.3;
1435       Double_t sy3=1.5;
1436       //sz3=1.5;
1437       //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
1438       //      Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
1439       //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
1440
1441       Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
1442       Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
1443       Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
1444       Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
1445       Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
1446       Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
1447       Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
1448       Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
1449       Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
1450       Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
1451
1452       c[0]=sy1;
1453       c[1]=0.;       c[2]=sz1;
1454       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
1455       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
1456       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
1457       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
1458       c[13]=f30*sy1*f40+f32*sy2*f42;
1459       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
1460
1461       UInt_t index=0;
1462       //kr0.GetIndex(is);
1463       AliTPCseed *track=new AliTPCseed(index, x, c, x1, sec*alpha+shift);
1464       track->fStopped =kFALSE;
1465       track->fIsSeeding = kTRUE;
1466       Int_t rc=FollowProlongation(*track, i2);
1467         track->fLastPoint = i1;  // first cluster in track position
1468         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
1469         else arr->AddLast(track);
1470     }
1471   }
1472
1473
1474
1475 }
1476
1477
1478
1479
1480
1481 //_____________________________________________________________________________
1482 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
1483   //-----------------------------------------------------------------
1484   // This function reades track seeds.
1485   //-----------------------------------------------------------------
1486   TDirectory *savedir=gDirectory;
1487
1488   TFile *in=(TFile*)inp;
1489   if (!in->IsOpen()) {
1490      cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
1491      return 1;
1492   }
1493
1494   in->cd();
1495   TTree *seedTree=(TTree*)in->Get("Seeds");
1496   if (!seedTree) {
1497      cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
1498      cerr<<"can't get a tree with track seeds !\n";
1499      return 2;
1500   }
1501   AliTPCtrack *seed=new AliTPCtrack;
1502   seedTree->SetBranchAddress("tracks",&seed);
1503
1504   if (fSeeds==0) fSeeds=new TObjArray(15000);
1505
1506   Int_t n=(Int_t)seedTree->GetEntries();
1507   for (Int_t i=0; i<n; i++) {
1508      seedTree->GetEvent(i);
1509      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
1510   }
1511
1512   delete seed;
1513   delete seedTree;
1514   savedir->cd();
1515   return 0;
1516 }
1517
1518 //_____________________________________________________________________________
1519 Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
1520   //-----------------------------------------------------------------
1521   // This is a track finder.
1522   //-----------------------------------------------------------------
1523   TDirectory *savedir=gDirectory;
1524
1525   if (inp) {
1526      TFile *in=(TFile*)inp;
1527      if (!in->IsOpen()) {
1528         cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
1529         return 1;
1530      }
1531   }
1532
1533   if (!out->IsOpen()) {
1534      cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
1535      return 2;
1536   }
1537
1538   out->cd();
1539
1540   char   tname[100];
1541   sprintf(tname,"TreeT_TPC_%d",fEventN);
1542   TTree tracktree(tname,"Tree with TPC tracks");
1543   TTree seedtree("Seeds","Seeds");
1544   AliTPCtrack *iotrack=0;
1545   AliTPCseed  *ioseed=0;
1546   tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
1547   TStopwatch timer;
1548
1549   printf("Loading clusters \n");
1550   LoadClusters();
1551   printf("Time for loading clusters: \t");timer.Print();timer.Start();
1552
1553   printf("Loading outer sectors\n");
1554   LoadOuterSectors();
1555   printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
1556
1557   printf("Loading inner sectors\n");
1558   LoadInnerSectors();
1559   printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
1560   fSectors = fOuterSec;
1561   fN=fkNOS;
1562
1563
1564   //find track seeds
1565   MakeSeedsAll();
1566   printf("Time for seeding: \t"); timer.Print();timer.Start();
1567   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
1568   Int_t nrows=nlow+nup;
1569
1570   Int_t gap=Int_t(0.3*nrows);
1571   Int_t i;
1572   //RemoveOverlap(fSeeds,0.6,2);
1573   Int_t nseed=fSeeds->GetEntriesFast();
1574   // outer sectors parallel tracking
1575   ParallelTracking(fSectors->GetNRows()-gap-1,0);
1576   //ParallelTracking(fSectors->GetNRows()-1,0);
1577   //RemoveOverlap(fSeeds, 0.6,3);
1578   //  ParallelTracking(49,0);
1579   printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
1580
1581   RemoveOverlap(fSeeds, 0.6,3);
1582   printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
1583   //parallel tracking
1584   fSectors = fInnerSec;
1585   fN=fkNIS;
1586   ParallelTracking(fSectors->GetNRows()-1,0);
1587   printf("Number of tracks after  inner tracking  %d\n",fNtracks);
1588   printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
1589   //
1590   RemoveOverlap(fSeeds,0.6,5,kTRUE);  // remove overlap -  shared points signed
1591   printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
1592   //
1593   //
1594
1595   ioseed  = (AliTPCseed*)(fSeeds->UncheckedAt(0));
1596   AliTPCseed * vseed = new AliTPCseed;
1597   vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
1598   vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
1599   vseed->fPoints->ExpandCreateFast(2);
1600
1601   TBranch * seedbranch =   seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
1602   //delete vseed;
1603   nseed=fSeeds->GetEntriesFast();
1604
1605   Int_t found = 0;
1606   for (i=0; i<nseed; i++) {
1607     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1608     if (!pt) continue;
1609     Int_t nc=t.GetNumberOfClusters();
1610     if (nc<20) continue;
1611     t.CookdEdx(0.02,0.6);
1612     CookLabel(pt,0.1); //For comparison only
1613     //   if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
1614     if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (t.fNFoundable>Int_t(0.3*nrows)))){
1615       iotrack=pt;
1616       tracktree.Fill();
1617      cerr<<found++<<'\r';
1618     }
1619     else
1620       if ( (pt->IsActive())) fNtracks--;
1621     pt->RebuildSeed();
1622     seedbranch->SetAddress(&pt);
1623
1624     seedtree.Fill();
1625     for (Int_t j=0;j<160;j++){
1626       delete pt->fPoints->RemoveAt(j);
1627     }
1628     delete pt->fPoints;
1629     pt->fPoints =0;
1630     delete fSeeds->RemoveAt(i);
1631   }
1632   printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
1633
1634   UnloadClusters();
1635   printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
1636
1637   tracktree.Write();
1638   seedtree.Write();
1639   cerr<<"Number of found tracks : "<<fNtracks<<"\t"<<found<<endl;
1640
1641   savedir->cd();
1642
1643   return 0;
1644 }
1645
1646
1647 void  AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
1648 {
1649   //
1650   // try to track in parralel
1651
1652   Int_t nseed=fSeeds->GetEntriesFast();
1653   //prepare seeds for tracking
1654   for (Int_t i=0; i<nseed; i++) {
1655     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1656     if (!pt) continue;
1657     if (!t.IsActive()) continue;
1658     // follow prolongation to the first layer
1659     FollowProlongation(t, rfirst+1);
1660   }
1661
1662
1663   //
1664   for (Int_t nr=rfirst; nr>=rlast; nr--){
1665     // make indexes with the cluster tracks for given
1666     //    for (Int_t i = 0;i<fN;i++)
1667     //  fSectors[i][nr].MakeClusterTracks();
1668
1669     // find nearest cluster
1670     for (Int_t i=0; i<nseed; i++) {
1671       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
1672       if (!pt) continue;
1673       if (!pt->IsActive()) continue;
1674       if (pt->fRelativeSector>17) {
1675         continue;
1676       }
1677       UpdateClusters(t,i,nr);
1678     }
1679     // prolonagate to the nearest cluster - if founded
1680     for (Int_t i=0; i<nseed; i++) {
1681       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
1682       if (!pt) continue;
1683       if (!pt->IsActive()) continue;
1684       if (pt->fRelativeSector>17) {
1685         continue;
1686       }
1687       FollowToNextCluster(i,nr);
1688     }
1689     //    for (Int_t i= 0;i<fN;i++)
1690     //  fSectors[i][nr].ClearClusterTracks();
1691   }
1692
1693 }
1694
1695 Float_t  AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
1696 {
1697   //
1698   //
1699   Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1700   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1701   Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
1702   Float_t angular  = seed->GetSnp();
1703   angular = angular*angular/(1-angular*angular);
1704   //  angular*=angular;
1705   //angular  = TMath::Sqrt(angular/(1-angular));
1706   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
1707   return res;
1708 }
1709 Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
1710 {
1711   //
1712   //
1713   Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
1714   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
1715   Float_t sres = fParam->GetZSigma();
1716   Float_t angular  = seed->GetTgl();
1717   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
1718   return res;
1719 }
1720
1721
1722
1723 //_________________________________________________________________________
1724 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1725   //--------------------------------------------------------------------
1726   //       Return pointer to a given cluster
1727   //--------------------------------------------------------------------
1728   Int_t sec=(index&0xff000000)>>24;
1729   Int_t row=(index&0x00ff0000)>>16;
1730   Int_t ncl=(index&0x0000ffff)>>00;
1731
1732   AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
1733   if (!clrow) return 0;
1734   return (AliTPCclusterMI*)(*clrow)[ncl];
1735 }
1736
1737 //__________________________________________________________________________
1738 void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
1739   //--------------------------------------------------------------------
1740   //This function "cooks" a track label. If label<0, this track is fake.
1741   //--------------------------------------------------------------------
1742   Int_t noc=t->GetNumberOfClusters();
1743   Int_t *lb=new Int_t[noc];
1744   Int_t *mx=new Int_t[noc];
1745   AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
1746
1747   Int_t i;
1748   for (i=0; i<noc; i++) {
1749      lb[i]=mx[i]=0;
1750      Int_t index=t->GetClusterIndex(i);
1751      clusters[i]=GetClusterMI(index);
1752   }
1753
1754   Int_t lab=123456789;
1755   for (i=0; i<noc; i++) {
1756     AliTPCclusterMI *c=clusters[i];
1757     if (!clusters[i]) continue;
1758     lab=TMath::Abs(c->GetLabel(0));
1759     Int_t j;
1760     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
1761     lb[j]=lab;
1762     (mx[j])++;
1763   }
1764
1765   Int_t max=0;
1766   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
1767
1768   for (i=0; i<noc; i++) {
1769     AliTPCclusterMI *c=clusters[i];
1770     if (!clusters[i]) continue;
1771     if (TMath::Abs(c->GetLabel(1)) == lab ||
1772         TMath::Abs(c->GetLabel(2)) == lab ) max++;
1773   }
1774
1775   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
1776
1777   else {
1778      Int_t tail=Int_t(0.10*noc);
1779      max=0;
1780      for (i=1; i<=tail; i++) {
1781        AliTPCclusterMI *c=clusters[noc-i];
1782        if (!clusters[i]) continue;
1783        if (lab == TMath::Abs(c->GetLabel(0)) ||
1784            lab == TMath::Abs(c->GetLabel(1)) ||
1785            lab == TMath::Abs(c->GetLabel(2))) max++;
1786      }
1787      if (max < Int_t(0.5*tail)) lab=-lab;
1788   }
1789
1790   t->SetLabel(lab);
1791
1792   delete[] lb;
1793   delete[] mx;
1794   delete[] clusters;
1795 }
1796
1797 //_________________________________________________________________________
1798 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
1799   //-----------------------------------------------------------------------
1800   // Setup inner sector
1801   //-----------------------------------------------------------------------
1802   if (f==0) {
1803      fAlpha=par->GetInnerAngle();
1804      fAlphaShift=par->GetInnerAngleShift();
1805      fPadPitchWidth=par->GetInnerPadPitchWidth();
1806      fPadPitchLength=par->GetInnerPadPitchLength();
1807      fN=par->GetNRowLow();
1808      fRow=new AliTPCRow[fN];
1809      for (Int_t i=0; i<fN; i++) {
1810        fRow[i].SetX(par->GetPadRowRadiiLow(i));
1811        fRow[i].fDeadZone =1.5;  //1.5 cm of dead zone
1812      }
1813   } else {
1814      fAlpha=par->GetOuterAngle();
1815      fAlphaShift=par->GetOuterAngleShift();
1816      fPadPitchWidth  = par->GetOuterPadPitchWidth();
1817      fPadPitchLength = par->GetOuter1PadPitchLength();
1818      f1PadPitchLength = par->GetOuter1PadPitchLength();
1819      f2PadPitchLength = par->GetOuter2PadPitchLength();
1820
1821      fN=par->GetNRowUp();
1822      fRow=new AliTPCRow[fN];
1823      for (Int_t i=0; i<fN; i++) {
1824        fRow[i].SetX(par->GetPadRowRadiiUp(i));
1825        fRow[i].fDeadZone =1.5;  // 1.5 cm of dead zone
1826      }
1827   }
1828 }
1829
1830
1831 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
1832   //
1833   if (fClusterTracks) delete [] fClusterTracks;
1834   fClusterTracks = 0;
1835 }
1836
1837 void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
1838   //create cluster tracks
1839   if (fN>0)
1840     fClusterTracks = new AliTPCclusterTracks[fN];
1841 }
1842
1843 void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
1844   if (fClusterTracks) delete[] fClusterTracks;
1845   fClusterTracks =0;
1846 }
1847
1848
1849
1850 void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
1851   //
1852   //
1853   // update information of the cluster tracks - if track is nearer then other tracks to the
1854   // given track
1855   const AliTPCclusterMI * cl = (*this)[clindex];
1856   AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
1857   // find the distance of the cluster to the track
1858   Float_t dy2 = (cl->GetY()- seed->GetY());
1859   dy2*=dy2;
1860   Float_t dz2 = (cl->GetZ()- seed->GetZ());
1861   dz2*=dz2;
1862   //
1863   Float_t distance = TMath::Sqrt(dy2+dz2);
1864   if (distance> 3)
1865     return;  // MI - to be changed - AliTPCtrackerParam
1866
1867   if ( distance < cltracks->fDistance[0]){
1868     cltracks->fDistance[2] =cltracks->fDistance[1];
1869     cltracks->fDistance[1] =cltracks->fDistance[0];
1870     cltracks->fDistance[0] =distance;
1871     cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1872     cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
1873     cltracks->fTrackIndex[0] =trindex;
1874   }
1875   else
1876     if ( distance < cltracks->fDistance[1]){
1877       cltracks->fDistance[2] =cltracks->fDistance[1];
1878       cltracks->fDistance[1] =distance;
1879       cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
1880       cltracks->fTrackIndex[1] =trindex;
1881     } else
1882       if (distance < cltracks->fDistance[2]){
1883         cltracks->fDistance[2] =distance;
1884         cltracks->fTrackIndex[2] =trindex;
1885       }
1886 }
1887
1888
1889 //_________________________________________________________________________
1890 void
1891 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
1892   //-----------------------------------------------------------------------
1893   // Insert a cluster into this pad row in accordence with its y-coordinate
1894   //-----------------------------------------------------------------------
1895   if (fN==kMaxClusterPerRow) {
1896     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
1897   }
1898   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
1899   Int_t i=Find(c->GetZ());
1900   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
1901   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
1902   fIndex[i]=index; fClusters[i]=c; fN++;
1903 }
1904
1905 //___________________________________________________________________
1906 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
1907   //-----------------------------------------------------------------------
1908   // Return the index of the nearest cluster
1909   //-----------------------------------------------------------------------
1910   if (fN==0) return 0;
1911   if (z <= fClusters[0]->GetZ()) return 0;
1912   if (z > fClusters[fN-1]->GetZ()) return fN;
1913   Int_t b=0, e=fN-1, m=(b+e)/2;
1914   for (; b<e; m=(b+e)/2) {
1915     if (z > fClusters[m]->GetZ()) b=m+1;
1916     else e=m;
1917   }
1918   return m;
1919 }
1920
1921
1922
1923 //___________________________________________________________________
1924 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
1925   //-----------------------------------------------------------------------
1926   // Return the index of the nearest cluster in z y
1927   //-----------------------------------------------------------------------
1928   Float_t maxdistance = roady*roady + roadz*roadz;
1929
1930   AliTPCclusterMI *cl =0;
1931   for (Int_t i=Find(z-roadz); i<fN; i++) {
1932       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
1933       if (c->GetZ() > z+roadz) break;
1934       if ( (c->GetY()-y) >  roady ) continue;
1935       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
1936       if (maxdistance>distance) {
1937         maxdistance = distance;
1938         cl=c;
1939       }
1940   }
1941   return cl;
1942 }
1943
1944
1945
1946
1947
1948 AliTPCseed::AliTPCseed():AliTPCtrack(){
1949   //
1950   fRow=0;
1951   fRemoval =0;
1952   memset(fClusterIndex,0,sizeof(Int_t)*200);
1953   fPoints = 0;
1954   fEPoints = 0;
1955   fNFoundable =0;
1956   fNShared  =0;
1957   fTrackPoints =0;
1958   fRemoval = 0;
1959 }
1960
1961 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
1962   fPoints = 0;
1963   fEPoints = 0;
1964   fNShared  =0;
1965   fTrackPoints =0;
1966   fRemoval =0;
1967 }
1968
1969 AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
1970   fRow=0;
1971   memset(fClusterIndex,0,sizeof(Int_t)*200);
1972   fPoints = 0;
1973   fEPoints = 0;
1974   fNFoundable =0;
1975   fNShared  =0;
1976   fTrackPoints =0;
1977   fRemoval =0;
1978 }
1979
1980 AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
1981                                         Double_t xr, Double_t alpha):
1982   AliTPCtrack(index, xx, cc, xr, alpha) {
1983   //
1984   //
1985   fRow =0;
1986   memset(fClusterIndex,0,sizeof(Int_t)*200);
1987   fPoints = 0;
1988   fEPoints = 0;
1989   fNFoundable =0;
1990   fNShared  = 0;
1991   fTrackPoints =0;
1992   fRemoval =0;
1993 }
1994
1995 AliTPCseed::~AliTPCseed(){
1996   if (fPoints) delete fPoints;
1997   fPoints =0;
1998   fEPoints = 0;
1999   if (fTrackPoints){
2000     for (Int_t i=0;i<8;i++){
2001       delete [] fTrackPoints[i];
2002     }
2003     delete fTrackPoints;
2004     fTrackPoints =0;
2005   }
2006
2007 }
2008
2009 AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
2010 {
2011   //
2012   //
2013   if (!fTrackPoints) {
2014     fTrackPoints = new AliTPCTrackPoint*[8];
2015     for ( Int_t i=0;i<8;i++)
2016       fTrackPoints[i]=0;
2017   }
2018   Int_t index1 = i/20;
2019   if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
2020   return &(fTrackPoints[index1][i%20]);
2021 }
2022
2023 void AliTPCseed::RebuildSeed()
2024 {
2025   //
2026   // rebuild seed to be ready for storing
2027   fPoints = new TClonesArray("AliTPCTrackPoint",160);
2028   fPoints->ExpandCreateFast(160);
2029   fEPoints = new TClonesArray("AliTPCExactPoint",1);
2030   for (Int_t i=0;i<160;i++){
2031     AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
2032     *trpoint = *(GetTrackPoint(i));
2033   }
2034
2035 }
2036
2037 //_____________________________________________________________________________
2038 void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
2039   //-----------------------------------------------------------------
2040   // This funtion calculates dE/dX within the "low" and "up" cuts.
2041   //-----------------------------------------------------------------
2042
2043   Float_t amp[200];
2044   Float_t angular[200];
2045   Float_t weight[200];
2046   Int_t index[200];
2047   //Int_t nc = 0;
2048   //  TClonesArray & arr = *fPoints;
2049   Float_t meanlog = 100.;
2050
2051   Float_t mean[4]  = {0,0,0,0};
2052   Float_t sigma[4] = {1000,1000,1000,1000};
2053   Int_t nc[4]      = {0,0,0,0};
2054   Float_t norm[4]    = {1000,1000,1000,1000};
2055   //
2056   //
2057   fNShared =0;
2058
2059   for (Int_t of =0; of<4; of++){
2060     for (Int_t i=of;i<160;i+=4)
2061       {
2062         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
2063         AliTPCTrackPoint * point = GetTrackPoint(i);
2064         if (point==0) continue;
2065         if (point->fIsShared){
2066           fNShared++;
2067           continue;
2068         }
2069         if (point->GetCPoint().GetMax()<5) continue;
2070         Float_t angley = point->GetTPoint().GetAngleY();
2071         Float_t anglez = point->GetTPoint().GetAngleZ();
2072         Int_t   type   = point->GetCPoint().GetType();
2073         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
2074         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
2075         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
2076
2077         Float_t ampc   = 0;     // normalization to the number of electrons
2078         if (i>64){
2079           ampc = 1.*point->GetCPoint().GetMax();
2080           //ampc = 1.*point->GetCPoint().GetQ();
2081           //      AliTPCClusterPoint & p = point->GetCPoint();
2082           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
2083           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2084           //Float_t dz =
2085           //  TMath::Abs( Int_t(iz) - iz + 0.5);
2086           //ampc *= 1.15*(1-0.3*dy);
2087           //ampc *= 1.15*(1-0.3*dz);
2088           //      Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
2089           //ampc               *=zfactor;
2090         }
2091         else{
2092           ampc = 1.0*point->GetCPoint().GetMax();
2093           //ampc = 1.0*point->GetCPoint().GetQ();
2094           //AliTPCClusterPoint & p = point->GetCPoint();
2095           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
2096           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
2097           //Float_t dz =
2098           //  TMath::Abs( Int_t(iz) - iz + 0.5);
2099
2100           //ampc *= 1.15*(1-0.3*dy);
2101           //ampc *= 1.15*(1-0.3*dz);
2102           //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
2103           //ampc               *=zfactor;
2104
2105         }
2106         ampc *= 2.0;     // put mean value to channel 50
2107         //ampc *= 0.58;     // put mean value to channel 50
2108         Float_t w      =  1.;
2109         //      if (type>0)  w =  1./(type/2.-0.5);
2110         Float_t z = TMath::Abs(point->GetCPoint().GetZ());
2111         if (i<64) {
2112           ampc /= 0.6;
2113           //ampc /= (1+0.0008*z);
2114         } else
2115           if (i>128){
2116             ampc /=1.5;
2117             //ampc /= (1+0.0008*z);
2118           }else{
2119             //ampc /= (1+0.0008*z);
2120           }
2121
2122         if (type<0) {  //amp at the border - lower weight
2123           // w*= 2.;
2124
2125           continue;
2126         }
2127         if (rsigma>1.5) ampc/=1.3;  // if big backround
2128         amp[nc[of]]        = ampc;
2129         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
2130         weight[nc[of]]     = w;
2131         nc[of]++;
2132       }
2133
2134     TMath::Sort(nc[of],amp,index,kFALSE);
2135     Float_t sumamp=0;
2136     Float_t sumamp2=0;
2137     Float_t sumw=0;
2138     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
2139     meanlog = 200;
2140     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
2141       Float_t ampl      = amp[index[i]]/angular[index[i]];
2142       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
2143       //
2144       sumw    += weight[index[i]];
2145       sumamp  += weight[index[i]]*ampl;
2146       sumamp2 += weight[index[i]]*ampl*ampl;
2147       norm[of]    += angular[index[i]]*weight[index[i]];
2148     }
2149     if (sumw<1){
2150       SetdEdx(0);
2151     }
2152     else {
2153       norm[of] /= sumw;
2154       mean[of]  = sumamp/sumw;
2155       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
2156       if (sigma[of]>0.1)
2157         sigma[of] = TMath::Sqrt(sigma[of]);
2158       else
2159         sigma[of] = 1000;
2160
2161     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
2162     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
2163     //mean *=(1-0.1*(norm-1.));
2164     }
2165   }
2166
2167   Float_t dedx =0;
2168   fSdEdx =0;
2169   fMAngular =0;
2170   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
2171   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
2172
2173
2174   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
2175   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
2176
2177   Int_t norm2 = 0;
2178   Int_t norm3 = 0;
2179   for (Int_t i =0;i<4;i++){
2180     if (nc[i]>2&&nc[i]<1000){
2181       dedx      += mean[i] *nc[i];
2182       fSdEdx    += sigma[i]*(nc[i]-2);
2183       fMAngular += norm[i] *nc[i];
2184       norm2     += nc[i];
2185       norm3     += nc[i]-2;
2186     }
2187     fDEDX[i]  = mean[i];
2188     fSDEDX[i] = sigma[i];
2189     fNCDEDX[i]= nc[i];
2190   }
2191
2192   if (norm3>0){
2193     dedx   /=norm2;
2194     fSdEdx /=norm3;
2195     fMAngular/=norm2;
2196   }
2197   else{
2198     SetdEdx(0);
2199     return;
2200   }
2201   //  Float_t dedx1 =dedx;
2202
2203   dedx =0;
2204   for (Int_t i =0;i<4;i++){
2205     if (nc[i]>2&&nc[i]<1000){
2206       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
2207       dedx      += mean[i] *nc[i];
2208     }
2209     fDEDX[i]  = mean[i];
2210   }
2211   dedx /= norm2;
2212
2213
2214
2215   SetdEdx(dedx);
2216
2217   //mi deDX
2218
2219
2220
2221   //Very rough PID
2222   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
2223
2224   if (p<0.6) {
2225     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2226     if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
2227     SetMass(0.93827); return;
2228   }
2229
2230   if (p<1.2) {
2231     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
2232     SetMass(0.93827); return;
2233   }
2234
2235   SetMass(0.13957); return;
2236
2237 }
2238
2239