TPC/AliTPCtrackerMI.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17 $Log$
  18 Revision 1.15  2003/09/29 11:56:58  kowal2
  19 bug fix2
  20
  21 Revision 1.14  2003/09/29 11:39:43  kowal2
  22 bug fix
  23
  24 Revision 1.13  2003/09/29 11:28:19  kowal2
  25 completly rewritten
  26
  27 Revision 1.9.4.3  2003/06/23 14:47:10  hristov
  28 Minor fix
  29
  30 Revision 1.9.4.2  2003/06/23 10:06:13  hristov
  31 Updated information about the overlapping clusters (M.Ivanov)
  32
  33 Revision 1.9.4.1  2003/06/19 06:59:58  hristov
  34 Updated version of parallel tracking (M.Ivanov)
  35
  36 Revision 1.9  2003/03/19 17:14:11  hristov
  37 Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
  38
  39 Revision 1.8  2003/03/05 11:16:15  kowal2
  40 Logs added
  41
  42 */
  43
  44
  45
  46
  47
  48
  49
  50 /*
  51   AliTPC parallel tracker -
  52   How to use?  -
  53   run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
  54   run AliTPCFindTracksMI.C macro - to find tracks
  55   tracks are written to AliTPCtracks.root file
  56   for comparison also seeds are written to the same file - to special branch
  57 */
  58
  59 //-------------------------------------------------------
  60 //          Implementation of the TPC tracker
  61 //
  62 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
  63 //
  64 //-------------------------------------------------------
  65 #include <TObjArray.h>
  66 #include <TFile.h>
  67 #include <TTree.h>
  68 #include "Riostream.h"
  69
  70 #include "AliTPCtrackerMI.h"
  71 #include "AliTPCclusterMI.h"
  72 #include "AliTPCParam.h"
  73 #include "AliTPCClustersRow.h"
  74 #include "AliComplexCluster.h"
  75 #include "AliTPCpolyTrack.h"
  76 #include "TStopwatch.h"
  77 #include "AliESD.h"
  78 #include "AliHelix.h"
  79 #include "TGraph.h"
  80 //
  81 #include "AliRunLoader.h"
  82
  83
  84
  85 ClassImp(AliTPCseed)
  86 ClassImp(AliTPCtrackerMI)
  87
  88
  89 class TPCFastMath {
  90 public:
  91   TPCFastMath();
  92   static Double_t FastAsin(Double_t x);
  93  private:
  94   static Double_t fgFastAsin[20000];
  95 };
  96
  97 Double_t TPCFastMath::fgFastAsin[20000];
  98
  99 TPCFastMath::TPCFastMath(){
 100   for (Int_t i=0;i<10000;i++){
 101     fgFastAsin[2*i] = TMath::ASin(i/10000.);
 102     fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
 103   }
 104 }
 105
 106 Double_t TPCFastMath::FastAsin(Double_t x){
 107   if (x>0){
 108     Int_t index = int(x*10000);
 109     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
 110   }
 111   x*=-1;
 112   Int_t index = int(x*10000);
 113   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
 114 }
 115
 116 TPCFastMath gTPCMath;
 117
 118
 119
 120 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
 121
 122   AliTPCclusterMI* c =track->fCurrentCluster;
 123   if (accept>0) track->fCurrentClusterIndex1 |=0x8000;  //sign not accepted clusters
 124
 125   UInt_t i = track->fCurrentClusterIndex1;
 126
 127   Int_t sec=(i&0xff000000)>>24;
 128   //Int_t row = (i&0x00ff0000)>>16;
 129   track->fRow=(i&0x00ff0000)>>16;
 130   track->fSector = sec;
 131   //  Int_t index = i&0xFFFF;
 132   if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow();
 133   track->SetClusterIndex2(track->fRow, i);
 134
 135   //track->fFirstPoint = row;
 136   //if ( track->fLastPoint<row) track->fLastPoint =row;
 137   //  if (track->fRow<0 || track->fRow>160) {
 138   //  printf("problem\n");
 139   //}
 140   if (track->fFirstPoint>track->fRow)
 141     track->fFirstPoint = track->fRow;
 142   if (track->fLastPoint<track->fRow)
 143     track->fLastPoint  = track->fRow;
 144
 145
 146   track->fClusterPointer[track->fRow] = c;
 147   //
 148
 149   Float_t angle2 = track->GetSnp()*track->GetSnp();
 150   angle2 = TMath::Sqrt(angle2/(1-angle2));
 151   //
 152   //SET NEW Track Point
 153   //
 154   //  if (debug)
 155   {
 156     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->fRow));
 157     //
 158     point.SetSigmaY(c->GetSigmaY2()/track->fCurrentSigmaY2);
 159     point.SetSigmaZ(c->GetSigmaZ2()/track->fCurrentSigmaZ2);
 160     point.SetErrY(sqrt(track->fErrorY2));
 161     point.SetErrZ(sqrt(track->fErrorZ2));
 162     //
 163     point.SetX(track->GetX());
 164     point.SetY(track->GetY());
 165     point.SetZ(track->GetZ());
 166     point.SetAngleY(angle2);
 167     point.SetAngleZ(track->GetTgl());
 168     if (point.fIsShared){
 169       track->fErrorY2 *= 4;
 170       track->fErrorZ2 *= 4;
 171     }
 172   }
 173
 174   Double_t chi2 = track->GetPredictedChi2(track->fCurrentCluster);
 175   //
 176   track->fErrorY2 *= 1.3;
 177   track->fErrorY2 += 0.01;
 178   track->fErrorZ2 *= 1.3;
 179   track->fErrorZ2 += 0.005;
 180     //}
 181   if (accept>0) return 0;
 182   if (track->GetNumberOfClusters()%20==0){
 183     //    if (track->fHelixIn){
 184     //  TClonesArray & larr = *(track->fHelixIn);
 185     //  Int_t ihelix = larr.GetEntriesFast();
 186     //  new(larr[ihelix]) AliHelix(*track) ;
 187     //}
 188   }
 189   track->fNoCluster =0;
 190   return track->Update(c,chi2,i);
 191 }
 192
 193
 194
 195 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster, Float_t factor,
 196                                       Float_t cory, Float_t corz)
 197 {
 198   //
 199   // decide according desired precision to accept given
 200   // cluster for tracking
 201   Double_t sy2=ErrY2(seed,cluster)*cory;
 202   Double_t sz2=ErrZ2(seed,cluster)*corz;
 203   //sy2=ErrY2(seed,cluster)*cory;
 204   //sz2=ErrZ2(seed,cluster)*cory;
 205
 206   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
 207   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
 208
 209   Double_t rdistancey2 = (seed->fCurrentCluster->GetY()-seed->GetY())*
 210     (seed->fCurrentCluster->GetY()-seed->GetY())/sdistancey2;
 211   Double_t rdistancez2 = (seed->fCurrentCluster->GetZ()-seed->GetZ())*
 212     (seed->fCurrentCluster->GetZ()-seed->GetZ())/sdistancez2;
 213
 214   Double_t rdistance2  = rdistancey2+rdistancez2;
 215   //Int_t  accept =0;
 216
 217   if (rdistance2>16) return 3;
 218
 219
 220   if ((rdistancey2>9.*factor || rdistancez2>9.*factor) && cluster->GetType()==0)
 221     return 2;  //suspisiouce - will be changed
 222
 223   if ((rdistancey2>6.25*factor || rdistancez2>6.25*factor) && cluster->GetType()>0)
 224     // strict cut on overlaped cluster
 225     return  2;  //suspisiouce - will be changed
 226
 227   if ( (rdistancey2>1.*factor || rdistancez2>6.25*factor )
 228        && cluster->GetType()<0){
 229     seed->fNFoundable--;
 230     return 2;
 231   }
 232   return 0;
 233 }
 234
 235
 236
 237
 238 //_____________________________________________________________________________
 239 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
 240 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 241 {
 242   //---------------------------------------------------------------------
 243   // The main TPC tracker constructor
 244   //---------------------------------------------------------------------
 245   fInnerSec=new AliTPCSector[fkNIS];
 246   fOuterSec=new AliTPCSector[fkNOS];
 247
 248   Int_t i;
 249   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
 250   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
 251
 252   fN=0;  fSectors=0;
 253
 254   fSeeds=0;
 255   fNtracks = 0;
 256   fParam = par;
 257   Int_t nrowlow = par->GetNRowLow();
 258   Int_t nrowup = par->GetNRowUp();
 259
 260
 261   for (Int_t i=0;i<nrowlow;i++){
 262     fXRow[i]     = par->GetPadRowRadiiLow(i);
 263     fPadLength[i]= par->GetPadPitchLength(0,i);
 264     fYMax[i]     = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
 265   }
 266
 267
 268   for (Int_t i=0;i<nrowup;i++){
 269     fXRow[i+nrowlow]      = par->GetPadRowRadiiUp(i);
 270     fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
 271     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
 272   }
 273   fSeeds=0;
 274   //
 275   fInput    = 0;
 276   fOutput   = 0;
 277   fSeedTree = 0;
 278   fTreeDebug =0;
 279   fNewIO     =0;
 280   fDebug     =0;
 281   fEvent     =0;
 282 }
 283
 284 //_____________________________________________________________________________
 285 AliTPCtrackerMI::~AliTPCtrackerMI() {
 286   //------------------------------------------------------------------
 287   // TPC tracker destructor
 288   //------------------------------------------------------------------
 289   delete[] fInnerSec;
 290   delete[] fOuterSec;
 291   if (fSeeds) {
 292     fSeeds->Delete();
 293     delete fSeeds;
 294   }
 295 }
 296
 297 void AliTPCtrackerMI::SetIO()
 298 {
 299   //
 300   fNewIO   =  kTRUE;
 301   fInput   =  AliRunLoader::GetTreeR("TPC", kFALSE,AliConfig::fgkDefaultEventFolderName);
 302   fOutput  =  AliRunLoader::GetTreeT("TPC", kTRUE,AliConfig::fgkDefaultEventFolderName);
 303   AliTPCtrack *iotrack= new AliTPCtrack;
 304   //    iotrack->fHelixIn   = new TClonesArray("AliHelix");
 305   //iotrack->fHelixOut  = new TClonesArray("AliHelix");
 306   fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 307   delete iotrack;
 308 }
 309
 310 void AliTPCtrackerMI::SetIO(TTree * input, TTree * output, AliESD * event)
 311 {
 312
 313   // set input
 314   fNewIO = kFALSE;
 315   fInput    = 0;
 316   fOutput   = 0;
 317   fSeedTree = 0;
 318   fTreeDebug =0;
 319   fInput = input;
 320   if (input==0){
 321     return;
 322   }
 323   //set output
 324   fOutput = output;
 325   if (output){
 326     AliTPCtrack *iotrack= new AliTPCtrack;
 327     //    iotrack->fHelixIn   = new TClonesArray("AliHelix");
 328     //iotrack->fHelixOut  = new TClonesArray("AliHelix");
 329     fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 330     delete iotrack;
 331   }
 332   if (output && (fDebug&2)){
 333     //write the full seed information if specified in debug mode
 334     //
 335     fSeedTree =  new TTree("Seeds","Seeds");
 336     AliTPCseed * vseed = new AliTPCseed;
 337     //
 338     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
 339     arrtr->ExpandCreateFast(160);
 340     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
 341     //
 342     vseed->fPoints = arrtr;
 343     vseed->fEPoints = arre;
 344     //    vseed->fClusterPoints = arrcl;
 345     fSeedTree->Branch("seeds","AliTPCseed",&vseed,32000,99);
 346     delete arrtr;
 347     delete arre;
 348     fTreeDebug = new TTree("trackDebug","trackDebug");
 349     TClonesArray * arrd = new TClonesArray("AliTPCTrackPoint2",0);
 350     fTreeDebug->Branch("debug",&arrd,32000,99);
 351   }
 352
 353
 354   //set ESD event
 355   fEvent  = event;
 356 }
 357
 358 void AliTPCtrackerMI::WriteTracks()
 359 {
 360   //
 361   // write tracks to the given output tree -
 362   // output specified with SetIO routine
 363   if (!fSeeds)  return;
 364   if (fEvent){
 365     // write tracks to the event
 366     // store index of the track
 367     Int_t nseed=fSeeds->GetEntriesFast();
 368     for (Int_t i=0; i<nseed; i++) {
 369       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
 370       if (!pt) continue;
 371       AliESDtrack iotrack;
 372       iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 373       //iotrack.SetTPCindex(i);
 374       fEvent->AddTrack(&iotrack);
 375     }
 376   }
 377
 378
 379   if (fOutput){
 380     AliTPCtrack *iotrack= 0;
 381     Int_t nseed=fSeeds->GetEntriesFast();
 382     for (Int_t i=0; i<nseed; i++) {
 383       iotrack= (AliTPCtrack*)fSeeds->UncheckedAt(i);
 384       if (iotrack) break;
 385     }
 386
 387     //TBranch * br = fOutput->Branch("tracks","AliTPCtrack",&iotrack,32000,100);
 388     TBranch * br = fOutput->GetBranch("tracks");
 389     br->SetAddress(&iotrack);
 390     //
 391     for (Int_t i=0; i<nseed; i++) {
 392       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
 393       if (!pt) continue;
 394       iotrack = pt;
 395       pt->fLab2 =i;
 396       //      br->SetAddress(&iotrack);
 397       fOutput->Fill();
 398       iotrack =0;
 399     }
 400   }
 401     // delete iotrack;
 402     //
 403   if (fSeedTree){
 404     //write the full seed information if specified in debug mode
 405
 406     AliTPCseed * vseed = new AliTPCseed;
 407     //
 408     TClonesArray * arrtr = new TClonesArray("AliTPCTrackPoint",160);
 409     arrtr->ExpandCreateFast(160);
 410     //TClonesArray * arrcl = new TClonesArray("AliTPCclusterMI",160);
 411     //arrcl->ExpandCreateFast(160);
 412     TClonesArray * arre = new TClonesArray("AliTPCExactPoint",160);
 413     //
 414     vseed->fPoints = arrtr;
 415     vseed->fEPoints = arre;
 416     //    vseed->fClusterPoints = arrcl;
 417     //TBranch * brseed = seedtree->Branch("seeds","AliTPCseed",&vseed,32000,99);
 418     TBranch * brseed = fSeedTree->GetBranch("seeds");
 419
 420     Int_t nseed=fSeeds->GetEntriesFast();
 421
 422     for (Int_t i=0; i<nseed; i++) {
 423       AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i);
 424       if (!pt) continue;
 425       pt->fPoints = arrtr;
 426       //      pt->fClusterPoints = arrcl;
 427       pt->fEPoints       = arre;
 428       pt->RebuildSeed();
 429       vseed = pt;
 430       brseed->SetAddress(&vseed);
 431       fSeedTree->Fill();
 432       pt->fPoints  = 0;
 433       pt->fEPoints = 0;
 434       //      pt->fClusterPoints = 0;
 435     }
 436     fSeedTree->Write();
 437     if (fTreeDebug) fTreeDebug->Write();
 438   }
 439
 440 }
 441
 442
 443
 444
 445 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
 446   //
 447   //
 448   //seed->SetErrorY2(0.1);
 449   //return 0.1;
 450   //calculate look-up table at the beginning
 451   static Bool_t  ginit = kFALSE;
 452   static Float_t gnoise1,gnoise2,gnoise3;
 453   static Float_t ggg1[10000];
 454   static Float_t ggg2[10000];
 455   static Float_t ggg3[10000];
 456   static Float_t glandau1[10000];
 457   static Float_t glandau2[10000];
 458   static Float_t glandau3[10000];
 459   //
 460   static Float_t gcor01[500];
 461   static Float_t gcor02[500];
 462   static Float_t gcorp[500];
 463   //
 464
 465   //
 466   if (ginit==kFALSE){
 467     for (Int_t i=1;i<500;i++){
 468       Float_t rsigma = float(i)/100.;
 469       gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
 470       gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
 471       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 472     }
 473
 474     //
 475     for (Int_t i=3;i<10000;i++){
 476       //
 477       //
 478       // inner sector
 479       Float_t amp = float(i);
 480       Float_t padlength =0.75;
 481       gnoise1 = 0.0004/padlength;
 482       Float_t nel     = 0.268*amp;
 483       Float_t nprim   = 0.155*amp;
 484       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 485       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 486       if (glandau1[i]>1) glandau1[i]=1;
 487       glandau1[i]*=padlength*padlength/12.;
 488       //
 489       // outer short
 490       padlength =1.;
 491       gnoise2   = 0.0004/padlength;
 492       nel       = 0.3*amp;
 493       nprim     = 0.133*amp;
 494       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 495       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 496       if (glandau2[i]>1) glandau2[i]=1;
 497       glandau2[i]*=padlength*padlength/12.;
 498       //
 499       //
 500       // outer long
 501       padlength =1.5;
 502       gnoise3   = 0.0004/padlength;
 503       nel       = 0.3*amp;
 504       nprim     = 0.133*amp;
 505       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 506       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 507       if (glandau3[i]>1) glandau3[i]=1;
 508       glandau3[i]*=padlength*padlength/12.;
 509       //
 510     }
 511     ginit = kTRUE;
 512   }
 513   //
 514   //
 515   //
 516   Int_t amp = int(TMath::Abs(cl->GetQ()));
 517   if (amp>9999) {
 518     seed->SetErrorY2(1.);
 519     return 1.;
 520   }
 521   Float_t snoise2;
 522   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 523   Int_t ctype = cl->GetType();
 524   Float_t padlength= GetPadPitchLength(seed->fRow);
 525   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 526   angle2 = angle2/(1-angle2);
 527   //
 528   //cluster "quality"
 529   Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->fCurrentSigmaY2));
 530   Float_t res;
 531   //
 532   if (fSectors==fInnerSec){
 533     snoise2 = gnoise1;
 534     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 535     if (ctype==0) res *= gcor01[rsigmay];
 536     if ((ctype>0)){
 537       res+=0.002;
 538       res*= gcorp[rsigmay];
 539     }
 540   }
 541   else {
 542     if (padlength<1.1){
 543       snoise2 = gnoise2;
 544       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 545       if (ctype==0) res *= gcor02[rsigmay];
 546       if ((ctype>0)){
 547         res+=0.002;
 548         res*= gcorp[rsigmay];
 549       }
 550     }
 551     else{
 552       snoise2 = gnoise3;
 553       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 554       if (ctype==0) res *= gcor02[rsigmay];
 555       if ((ctype>0)){
 556         res+=0.002;
 557         res*= gcorp[rsigmay];
 558       }
 559     }
 560   }
 561
 562   if (ctype<0){
 563     res+=0.005;
 564     res*=2.4;  // overestimate error 2 times
 565   }
 566   res+= snoise2;
 567
 568   if (res<2*snoise2)
 569     res = 2*snoise2;
 570
 571   seed->SetErrorY2(res);
 572   return res;
 573
 574
 575 }
 576
 577
 578
 579 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 580   //
 581   //
 582   //seed->SetErrorY2(0.1);
 583   //return 0.1;
 584   //calculate look-up table at the beginning
 585   static Bool_t  ginit = kFALSE;
 586   static Float_t gnoise1,gnoise2,gnoise3;
 587   static Float_t ggg1[10000];
 588   static Float_t ggg2[10000];
 589   static Float_t ggg3[10000];
 590   static Float_t glandau1[10000];
 591   static Float_t glandau2[10000];
 592   static Float_t glandau3[10000];
 593   //
 594   static Float_t gcor01[1000];
 595   static Float_t gcor02[1000];
 596   static Float_t gcorp[1000];
 597   //
 598
 599   //
 600   if (ginit==kFALSE){
 601     for (Int_t i=1;i<1000;i++){
 602       Float_t rsigma = float(i)/100.;
 603       gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
 604       gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
 605       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 606     }
 607
 608     //
 609     for (Int_t i=3;i<10000;i++){
 610       //
 611       //
 612       // inner sector
 613       Float_t amp = float(i);
 614       Float_t padlength =0.75;
 615       gnoise1 = 0.0004/padlength;
 616       Float_t nel     = 0.268*amp;
 617       Float_t nprim   = 0.155*amp;
 618       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 619       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 620       if (glandau1[i]>1) glandau1[i]=1;
 621       glandau1[i]*=padlength*padlength/12.;
 622       //
 623       // outer short
 624       padlength =1.;
 625       gnoise2   = 0.0004/padlength;
 626       nel       = 0.3*amp;
 627       nprim     = 0.133*amp;
 628       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 629       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 630       if (glandau2[i]>1) glandau2[i]=1;
 631       glandau2[i]*=padlength*padlength/12.;
 632       //
 633       //
 634       // outer long
 635       padlength =1.5;
 636       gnoise3   = 0.0004/padlength;
 637       nel       = 0.3*amp;
 638       nprim     = 0.133*amp;
 639       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 640       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 641       if (glandau3[i]>1) glandau3[i]=1;
 642       glandau3[i]*=padlength*padlength/12.;
 643       //
 644     }
 645     ginit = kTRUE;
 646   }
 647   //
 648   //
 649   //
 650   Int_t amp = int(TMath::Abs(cl->GetQ()));
 651   if (amp>9999) {
 652     seed->SetErrorY2(1.);
 653     return 1.;
 654   }
 655   Float_t snoise2;
 656   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 657   Int_t ctype = cl->GetType();
 658   Float_t padlength= GetPadPitchLength(seed->fRow);
 659   //
 660   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 661   //  if (angle2<0.6) angle2 = 0.6;
 662   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
 663   //
 664   //cluster "quality"
 665   Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2));
 666   Float_t res;
 667   //
 668   if (fSectors==fInnerSec){
 669     snoise2 = gnoise1;
 670     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 671     if (ctype==0) res *= gcor01[rsigmaz];
 672     if ((ctype>0)){
 673       res+=0.002;
 674       res*= gcorp[rsigmaz];
 675     }
 676   }
 677   else {
 678     if (padlength<1.1){
 679       snoise2 = gnoise2;
 680       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 681       if (ctype==0) res *= gcor02[rsigmaz];
 682       if ((ctype>0)){
 683         res+=0.002;
 684         res*= gcorp[rsigmaz];
 685       }
 686     }
 687     else{
 688       snoise2 = gnoise3;
 689       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 690       if (ctype==0) res *= gcor02[rsigmaz];
 691       if ((ctype>0)){
 692         res+=0.002;
 693         res*= gcorp[rsigmaz];
 694       }
 695     }
 696   }
 697
 698   if (ctype<0){
 699     res+=0.002;
 700     res*=1.3;
 701   }
 702   if ((ctype<0) &&amp<70){
 703     res+=0.002;
 704     res*=1.3;
 705   }
 706   res += snoise2;
 707   if (res<2*snoise2)
 708      res = 2*snoise2;
 709   if (res>3) res =3;
 710   seed->SetErrorZ2(res);
 711   return res;
 712 }
 713
 714
 715
 716 /*
 717 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 718   //
 719   //
 720   //seed->SetErrorZ2(0.1);
 721   //return 0.1;
 722
 723   Float_t snoise2;
 724   Float_t z = TMath::Abs(fParam->GetZLength()-TMath::Abs(seed->GetZ()));
 725   //
 726   Float_t rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ2);
 727   Int_t ctype = cl->GetType();
 728   Float_t amp = TMath::Abs(cl->GetQ());
 729
 730   Float_t nel;
 731   Float_t nprim;
 732   //
 733   Float_t landau=2 ;    //landau fluctuation part
 734   Float_t gg=2;         // gg fluctuation part
 735   Float_t padlength= GetPadPitchLength(seed->GetX());
 736
 737   if (fSectors==fInnerSec){
 738     snoise2 = 0.0004/padlength;
 739     nel     = 0.268*amp;
 740     nprim   = 0.155*amp;
 741     gg      = (2+0.001*nel/(padlength*padlength))/nel;
 742     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 743     if (landau>1) landau=1;
 744   }
 745   else {
 746     snoise2 = 0.0004/padlength;
 747     nel     = 0.3*amp;
 748     nprim   = 0.133*amp;
 749     gg      = (2+0.0008*nel/(padlength*padlength))/nel;
 750     landau  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 751     if (landau>1) landau=1;
 752   }
 753   Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
 754
 755   //
 756   Float_t angle2 = seed->GetSnp()*seed->GetSnp();
 757   angle2 = TMath::Sqrt((1-angle2));
 758   if (angle2<0.6) angle2 = 0.6;
 759   //angle2 = 1;
 760
 761   Float_t angle = seed->GetTgl()/angle2;
 762   Float_t angular = landau*angle*angle*padlength*padlength/12.;
 763   Float_t res = sdiff + angular;
 764
 765
 766   if ((ctype==0) && (fSectors ==fOuterSec))
 767     res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
 768
 769   if ((ctype==0) && (fSectors ==fInnerSec))
 770     res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
 771
 772   if ((ctype>0)){
 773     res+=0.005;
 774     res*= TMath::Power(rsigmaz+0.5,1.5);  //0.31+0.147*ctype;
 775   }
 776   if (ctype<0){
 777     res+=0.002;
 778     res*=1.3;
 779   }
 780   if ((ctype<0) &&amp<70){
 781     res+=0.002;
 782     res*=1.3;
 783   }
 784   res += snoise2;
 785   if (res<2*snoise2)
 786      res = 2*snoise2;
 787
 788   seed->SetErrorZ2(res);
 789   return res;
 790 }
 791 */
 792
 793
 794
 795 void AliTPCseed::Reset(Bool_t all)
 796 {
 797   //
 798   //
 799   SetNumberOfClusters(0);
 800   fNFoundable = 0;
 801   SetChi2(0);
 802   ResetCovariance();
 803   /*
 804   if (fTrackPoints){
 805     for (Int_t i=0;i<8;i++){
 806       delete [] fTrackPoints[i];
 807     }
 808     delete fTrackPoints;
 809     fTrackPoints =0;
 810   }
 811   */
 812
 813   if (all){
 814     for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
 815     for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
 816   }
 817
 818 }
 819
 820
 821 void AliTPCseed::Modify(Double_t factor)
 822 {
 823
 824   //------------------------------------------------------------------
 825   //This function makes a track forget its history :)
 826   //------------------------------------------------------------------
 827   if (factor<=0) {
 828     ResetCovariance();
 829     return;
 830   }
 831   fC00*=factor;
 832   fC10*=factor;  fC11*=factor;
 833   fC20*=factor;  fC21*=factor;  fC22*=factor;
 834   fC30*=factor;  fC31*=factor;  fC32*=factor;  fC33*=factor;
 835   fC40*=factor;  fC41*=factor;  fC42*=factor;  fC43*=factor;  fC44*=factor;
 836   SetNumberOfClusters(0);
 837   fNFoundable =0;
 838   SetChi2(0);
 839   fRemoval = 0;
 840   fCurrentSigmaY2 = 0.000005;
 841   fCurrentSigmaZ2 = 0.000005;
 842   fNoCluster     = 0;
 843   //fFirstPoint = 160;
 844   //fLastPoint  = 0;
 845 }
 846
 847
 848
 849
 850 Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
 851 {
 852   //-----------------------------------------------------------------
 853   // This function find proloncation of a track to a reference plane x=xk.
 854   // doesn't change internal state of the track
 855   //-----------------------------------------------------------------
 856
 857   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
 858
 859   if (TMath::Abs(fP4*xk - fP2) >= 0.999) {
 860     return 0;
 861   }
 862
 863   //  Double_t y1=fP0, z1=fP1;
 864   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
 865   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
 866
 867   y = fP0;
 868   z = fP1;
 869   //y += dx*(c1+c2)/(r1+r2);
 870   //z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
 871
 872   Double_t dy = dx*(c1+c2)/(r1+r2);
 873   Double_t dz = 0;
 874   //
 875   Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
 876   /*
 877   if (TMath::Abs(delta)>0.0001){
 878     dz = fP3*TMath::ASin(delta)/fP4;
 879   }else{
 880     dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
 881   }
 882   */
 883   dz =  fP3*TPCFastMath::FastAsin(delta)/fP4;
 884   //
 885   y+=dy;
 886   z+=dz;
 887
 888
 889   return 1;
 890 }
 891
 892
 893 //_____________________________________________________________________________
 894 Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const
 895 {
 896   //-----------------------------------------------------------------
 897   // This function calculates a predicted chi2 increment.
 898   //-----------------------------------------------------------------
 899   //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
 900   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
 901   r00+=fC00; r01+=fC10; r11+=fC11;
 902
 903   Double_t det=r00*r11 - r01*r01;
 904   if (TMath::Abs(det) < 1.e-10) {
 905     Int_t n=GetNumberOfClusters();
 906     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
 907     return 1e10;
 908   }
 909   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 910
 911   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
 912
 913   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
 914 }
 915
 916
 917 //_________________________________________________________________________________________
 918
 919
 920 Int_t AliTPCseed::Compare(const TObject *o) const {
 921   //-----------------------------------------------------------------
 922   // This function compares tracks according to the sector - for given sector according z
 923   //-----------------------------------------------------------------
 924   AliTPCseed *t=(AliTPCseed*)o;
 925
 926   if (fSort == 0){
 927     if (t->fRelativeSector>fRelativeSector) return -1;
 928     if (t->fRelativeSector<fRelativeSector) return 1;
 929     Double_t z2 = t->GetZ();
 930     Double_t z1 = GetZ();
 931     if (z2>z1) return 1;
 932     if (z2<z1) return -1;
 933     return 0;
 934   }
 935   else {
 936     Float_t f2 =1;
 937     f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
 938     if (t->fBConstrain) f2=1.2;
 939
 940     Float_t f1 =1;
 941     f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
 942
 943     if (fBConstrain)   f1=1.2;
 944
 945     if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
 946     else return +1;
 947   }
 948 }
 949
 950 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
 951 {
 952   //rotate to track "local coordinata
 953   Float_t x = seed->GetX();
 954   Float_t y = seed->GetY();
 955   Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
 956
 957   if (y > ymax) {
 958     seed->fRelativeSector= (seed->fRelativeSector+1) % fN;
 959     if (!seed->Rotate(fSectors->GetAlpha()))
 960       return;
 961   } else if (y <-ymax) {
 962     seed->fRelativeSector= (seed->fRelativeSector-1+fN) % fN;
 963     if (!seed->Rotate(-fSectors->GetAlpha()))
 964       return;
 965   }
 966
 967 }
 968
 969
 970
 971
 972 //_____________________________________________________________________________
 973 Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
 974   //-----------------------------------------------------------------
 975   // This function associates a cluster with this track.
 976   //-----------------------------------------------------------------
 977   Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
 978
 979   r00+=fC00; r01+=fC10; r11+=fC11;
 980   Double_t det=r00*r11 - r01*r01;
 981   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 982
 983   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
 984   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
 985   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
 986   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
 987   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
 988
 989   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
 990   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
 991   if (TMath::Abs(cur*fX-eta) >= 0.9) {
 992     return 0;
 993   }
 994
 995   fP0 += k00*dy + k01*dz;
 996   fP1 += k10*dy + k11*dz;
 997   fP2  = eta;
 998   fP3 += k30*dy + k31*dz;
 999   fP4  = cur;
1000
1001   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
1002   Double_t c12=fC21, c13=fC31, c14=fC41;
1003
1004   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
1005   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
1006   fC40-=k00*c04+k01*c14;
1007
1008   fC11-=k10*c01+k11*fC11;
1009   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
1010   fC41-=k10*c04+k11*c14;
1011
1012   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
1013   fC42-=k20*c04+k21*c14;
1014
1015   fC33-=k30*c03+k31*c13;
1016   fC43-=k40*c03+k41*c13;
1017
1018   fC44-=k40*c04+k41*c14;
1019
1020   Int_t n=GetNumberOfClusters();
1021   //  fIndex[n]=index;
1022   SetNumberOfClusters(n+1);
1023   SetChi2(GetChi2()+chisq);
1024
1025   return 1;
1026 }
1027
1028
1029
1030 //_____________________________________________________________________________
1031 Double_t AliTPCtrackerMI::f1old(Double_t x1,Double_t y1,
1032                    Double_t x2,Double_t y2,
1033                    Double_t x3,Double_t y3)
1034 {
1035   //-----------------------------------------------------------------
1036   // Initial approximation of the track curvature
1037   //-----------------------------------------------------------------
1038   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1039   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1040                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1041   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1042                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1043
1044   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1045   if ( xr*xr+yr*yr<=0.00000000000001) return 100;
1046   return -xr*yr/sqrt(xr*xr+yr*yr);
1047 }
1048
1049
1050
1051 //_____________________________________________________________________________
1052 Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
1053                    Double_t x2,Double_t y2,
1054                    Double_t x3,Double_t y3)
1055 {
1056   //-----------------------------------------------------------------
1057   // Initial approximation of the track curvature
1058   //-----------------------------------------------------------------
1059   x3 -=x1;
1060   x2 -=x1;
1061   y3 -=y1;
1062   y2 -=y1;
1063   //
1064   Double_t det = x3*y2-x2*y3;
1065   if (det==0) {
1066     return 100;
1067   }
1068   //
1069   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1070   Double_t x0 = x3*0.5-y3*u;
1071   Double_t y0 = y3*0.5+x3*u;
1072   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1073   if (det<0) c2*=-1;
1074   return c2;
1075 }
1076
1077
1078 Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
1079                    Double_t x2,Double_t y2,
1080                    Double_t x3,Double_t y3)
1081 {
1082   //-----------------------------------------------------------------
1083   // Initial approximation of the track curvature
1084   //-----------------------------------------------------------------
1085   x3 -=x1;
1086   x2 -=x1;
1087   y3 -=y1;
1088   y2 -=y1;
1089   //
1090   Double_t det = x3*y2-x2*y3;
1091   if (det==0) {
1092     return 100;
1093   }
1094   //
1095   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
1096   Double_t x0 = x3*0.5-y3*u;
1097   Double_t y0 = y3*0.5+x3*u;
1098   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
1099   if (det<0) c2*=-1;
1100   x0+=x1;
1101   x0*=c2;
1102   return x0;
1103 }
1104
1105
1106
1107 //_____________________________________________________________________________
1108 Double_t AliTPCtrackerMI::f2old(Double_t x1,Double_t y1,
1109                    Double_t x2,Double_t y2,
1110                    Double_t x3,Double_t y3)
1111 {
1112   //-----------------------------------------------------------------
1113   // Initial approximation of the track curvature times center of curvature
1114   //-----------------------------------------------------------------
1115   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1116   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1117                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1118   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1119                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1120
1121   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1122
1123   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1124 }
1125
1126 //_____________________________________________________________________________
1127 Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1,
1128                    Double_t x2,Double_t y2,
1129                    Double_t z1,Double_t z2)
1130 {
1131   //-----------------------------------------------------------------
1132   // Initial approximation of the tangent of the track dip angle
1133   //-----------------------------------------------------------------
1134   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1135 }
1136
1137
1138 Double_t AliTPCtrackerMI::f3n(Double_t x1,Double_t y1,
1139                    Double_t x2,Double_t y2,
1140                    Double_t z1,Double_t z2, Double_t c)
1141 {
1142   //-----------------------------------------------------------------
1143   // Initial approximation of the tangent of the track dip angle
1144   //-----------------------------------------------------------------
1145
1146   //  Double_t angle1;
1147
1148   //angle1    =  (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1149   //
1150   Double_t d  =  TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1151   if (TMath::Abs(d*c*0.5)>1) return 0;
1152   //  Double_t   angle2    =  TMath::ASin(d*c*0.5);
1153   //  Double_t   angle2    =  TPCFastMath::FastAsin(d*c*0.5);
1154   Double_t   angle2    = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): TPCFastMath::FastAsin(d*c*0.5);
1155
1156   angle2  = (z1-z2)*c/(angle2*2.);
1157   return angle2;
1158 }
1159
1160 Bool_t   AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1161 {//-----------------------------------------------------------------
1162   // This function find proloncation of a track to a reference plane x=x2.
1163   //-----------------------------------------------------------------
1164
1165   Double_t dx=x2-x1;
1166
1167   if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1168     return kFALSE;
1169   }
1170
1171   Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1172   Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1173   y = x[0];
1174   z = x[1];
1175
1176   Double_t dy = dx*(c1+c2)/(r1+r2);
1177   Double_t dz = 0;
1178   //
1179   Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1180
1181   if (TMath::Abs(delta)>0.01){
1182     dz = x[3]*TMath::ASin(delta)/x[4];
1183   }else{
1184     dz = x[3]*TPCFastMath::FastAsin(delta)/x[4];
1185   }
1186
1187   //dz = x[3]*TPCFastMath::FastAsin(delta)/x[4];
1188
1189   y+=dy;
1190   z+=dz;
1191
1192   return kTRUE;
1193 }
1194
1195
1196
1197 Int_t  AliTPCtrackerMI::LoadClusters()
1198 {
1199   //
1200   // load clusters to the memory
1201   AliTPCClustersRow *clrow= new AliTPCClustersRow;
1202   clrow->SetClass("AliTPCclusterMI");
1203   clrow->SetArray(0);
1204   clrow->GetArray()->ExpandCreateFast(10000);
1205   //
1206   //  TTree * tree = fClustersArray.GetTree();
1207
1208   TTree * tree = fInput;
1209   TBranch * br = tree->GetBranch("Segment");
1210   br->SetAddress(&clrow);
1211   //
1212   Int_t j=Int_t(tree->GetEntries());
1213   for (Int_t i=0; i<j; i++) {
1214     br->GetEntry(i);
1215     //
1216     Int_t sec,row;
1217     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1218     //
1219     AliTPCRow * tpcrow=0;
1220     Int_t left=0;
1221     if (sec<fkNIS*2){
1222       tpcrow = &(fInnerSec[sec%fkNIS][row]);
1223       left = sec/fkNIS;
1224     }
1225     else{
1226       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1227       left = (sec-fkNIS*2)/fkNOS;
1228     }
1229     if (left ==0){
1230       tpcrow->fN1 = clrow->GetArray()->GetEntriesFast();
1231       tpcrow->fClusters1 = new AliTPCclusterMI[tpcrow->fN1];
1232       for (Int_t i=0;i<tpcrow->fN1;i++)
1233         tpcrow->fClusters1[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1234     }
1235     if (left ==1){
1236       tpcrow->fN2 = clrow->GetArray()->GetEntriesFast();
1237       tpcrow->fClusters2 = new AliTPCclusterMI[tpcrow->fN2];
1238       for (Int_t i=0;i<tpcrow->fN2;i++)
1239         tpcrow->fClusters2[i] = *(AliTPCclusterMI*)(clrow->GetArray()->At(i));
1240     }
1241   }
1242   //
1243   delete clrow;
1244   LoadOuterSectors();
1245   LoadInnerSectors();
1246   return 0;
1247 }
1248
1249
1250 void AliTPCtrackerMI::UnloadClusters()
1251 {
1252   //
1253   // unload clusters from the memory
1254   //
1255   Int_t nrows = fOuterSec->GetNRows();
1256   for (Int_t sec = 0;sec<fkNOS;sec++)
1257     for (Int_t row = 0;row<nrows;row++){
1258       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1259       if (tpcrow){
1260         if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1261         if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1262       }
1263     }
1264   //
1265   nrows = fInnerSec->GetNRows();
1266   for (Int_t sec = 0;sec<fkNIS;sec++)
1267     for (Int_t row = 0;row<nrows;row++){
1268       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1269       if (tpcrow){
1270         if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1271         if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1272       }
1273     }
1274
1275   return ;
1276 }
1277
1278
1279 //_____________________________________________________________________________
1280 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1281   //-----------------------------------------------------------------
1282   // This function fills outer TPC sectors with clusters.
1283   //-----------------------------------------------------------------
1284   Int_t nrows = fOuterSec->GetNRows();
1285   UInt_t index=0;
1286   for (Int_t sec = 0;sec<fkNOS;sec++)
1287     for (Int_t row = 0;row<nrows;row++){
1288       AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1289       Int_t sec2 = sec+2*fkNIS;
1290       //left
1291       Int_t ncl = tpcrow->fN1;
1292       while (ncl--) {
1293         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1294         index=(((sec2<<8)+row)<<16)+ncl;
1295         tpcrow->InsertCluster(c,index);
1296       }
1297       //right
1298       ncl = tpcrow->fN2;
1299       while (ncl--) {
1300         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1301         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1302         tpcrow->InsertCluster(c,index);
1303       }
1304       //
1305       // write indexes for fast acces
1306       //
1307       for (Int_t i=0;i<510;i++)
1308         tpcrow->fFastCluster[i]=-1;
1309       for (Int_t i=0;i<tpcrow->GetN();i++){
1310         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1311         tpcrow->fFastCluster[zi]=i;  // write index
1312       }
1313       Int_t last = 0;
1314       for (Int_t i=0;i<510;i++){
1315         if (tpcrow->fFastCluster[i]<0)
1316           tpcrow->fFastCluster[i] = last;
1317         else
1318           last = tpcrow->fFastCluster[i];
1319       }
1320     }
1321   fN=fkNOS;
1322   fSectors=fOuterSec;
1323   return 0;
1324 }
1325
1326
1327 //_____________________________________________________________________________
1328 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1329   //-----------------------------------------------------------------
1330   // This function fills inner TPC sectors with clusters.
1331   //-----------------------------------------------------------------
1332   Int_t nrows = fInnerSec->GetNRows();
1333   UInt_t index=0;
1334   for (Int_t sec = 0;sec<fkNIS;sec++)
1335     for (Int_t row = 0;row<nrows;row++){
1336       AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1337       //
1338       //left
1339       Int_t ncl = tpcrow->fN1;
1340       while (ncl--) {
1341         AliTPCclusterMI *c= &(tpcrow->fClusters1[ncl]);
1342         index=(((sec<<8)+row)<<16)+ncl;
1343         tpcrow->InsertCluster(c,index);
1344       }
1345       //right
1346       ncl = tpcrow->fN2;
1347       while (ncl--) {
1348         AliTPCclusterMI *c= &(tpcrow->fClusters2[ncl]);
1349         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1350         tpcrow->InsertCluster(c,index);
1351       }
1352       //
1353       // write indexes for fast acces
1354       //
1355       for (Int_t i=0;i<510;i++)
1356         tpcrow->fFastCluster[i]=-1;
1357       for (Int_t i=0;i<tpcrow->GetN();i++){
1358         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1359         tpcrow->fFastCluster[zi]=i;  // write index
1360       }
1361       Int_t last = 0;
1362       for (Int_t i=0;i<510;i++){
1363         if (tpcrow->fFastCluster[i]<0)
1364           tpcrow->fFastCluster[i] = last;
1365         else
1366           last = tpcrow->fFastCluster[i];
1367       }
1368
1369     }
1370
1371   fN=fkNIS;
1372   fSectors=fInnerSec;
1373   return 0;
1374 }
1375
1376
1377
1378 //_________________________________________________________________________
1379 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1380   //--------------------------------------------------------------------
1381   //       Return pointer to a given cluster
1382   //--------------------------------------------------------------------
1383   Int_t sec=(index&0xff000000)>>24;
1384   Int_t row=(index&0x00ff0000)>>16;
1385   Int_t ncl=(index&0x0000ffff)>>00;
1386
1387   const AliTPCRow * tpcrow=0;
1388   AliTPCclusterMI * clrow =0;
1389   if (sec<fkNIS*2){
1390     tpcrow = &(fInnerSec[sec%fkNIS][row]);
1391     if (sec<fkNIS)
1392       clrow = tpcrow->fClusters1;
1393     else
1394       clrow = tpcrow->fClusters2;
1395   }
1396   else{
1397     tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1398     if (sec-2*fkNIS<fkNOS)
1399       clrow = tpcrow->fClusters1;
1400     else
1401       clrow = tpcrow->fClusters2;
1402   }
1403   if (tpcrow==0) return 0;
1404   if (tpcrow->GetN()<=ncl) return 0;
1405   //  return (AliTPCclusterMI*)(*tpcrow)[ncl];
1406   return &(clrow[ncl]);
1407
1408 }
1409
1410
1411
1412 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1413   //-----------------------------------------------------------------
1414   // This function tries to find a track prolongation to next pad row
1415   //-----------------------------------------------------------------
1416   //
1417   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1418
1419   //  if (t.GetRadius()>x+10 ) return 0;
1420   //  t.PropagateTo(x+0.02);
1421   //t.PropagateTo(x+0.01);
1422   if (!t.PropagateTo(x)) {
1423     t.fRemoval = 10;
1424     return 0;
1425   }
1426   //
1427   Double_t  y=t.GetY(), z=t.GetZ();
1428   if (TMath::Abs(y)>ymax){
1429     if (y > ymax) {
1430       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1431       if (!t.Rotate(fSectors->GetAlpha()))
1432         return 0;
1433     } else if (y <-ymax) {
1434       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1435       if (!t.Rotate(-fSectors->GetAlpha()))
1436         return 0;
1437     }
1438     if (!t.PropagateTo(x)) {
1439       return 0;
1440     }
1441     y=t.GetY();
1442   }
1443   //
1444   // update current shape info every 3 pad-row
1445   if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1446     //t.fCurrentSigmaY = GetSigmaY(&t);
1447     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1448     GetShape(&t,nr);
1449   }
1450   //
1451   AliTPCclusterMI *cl=0;
1452   UInt_t index=0;
1453
1454
1455   //Int_t nr2 = nr;
1456   if (t.GetClusterIndex2(nr)>0){
1457     //
1458     //cl = GetClusterMI(t.GetClusterIndex2(nr));
1459     index = t.GetClusterIndex2(nr);
1460     cl = t.fClusterPointer[nr];
1461     if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1462     t.fCurrentClusterIndex1 = index;
1463   }
1464
1465   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1466   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1467   Double_t  roady  =1.;
1468   Double_t  roadz = 1.;
1469   //
1470   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1471     t.fInDead = kTRUE;
1472     t.SetClusterIndex2(nr,-1);
1473     return 0;
1474   }
1475   else
1476     {
1477       if (TMath::Abs(z)<(1.05*x+10)) t.fNFoundable++;
1478       else
1479         return 0;
1480     }
1481   //calculate
1482   if (cl){
1483     t.fCurrentCluster = cl;
1484     t.fRow = nr;
1485     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1486     if (accept<3) {
1487       //if founded cluster is acceptible
1488       UpdateTrack(&t,accept);
1489       return 1;
1490     }
1491   }
1492
1493   if (krow) {
1494     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1495     cl = krow.FindNearest2(y,z,roady,roadz,index);
1496     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1497   }
1498   //  t.fNoCluster++;
1499
1500   if (cl) {
1501     t.fCurrentCluster = cl;
1502     t.fRow = nr;
1503     Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1504
1505     if (t.fCurrentCluster->IsUsed(10)){
1506       //
1507       //
1508
1509       t.fNShared++;
1510       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1511         t.fRemoval =10;
1512         return 0;
1513       }
1514     }
1515
1516     if (accept<3) UpdateTrack(&t,accept);
1517
1518   } else {
1519     if (t.fNFoundable*0.5 > t.GetNumberOfClusters()) t.fRemoval=10;
1520
1521   }
1522   return 1;
1523 }
1524
1525 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1526   //-----------------------------------------------------------------
1527   // This function tries to find a track prolongation to next pad row
1528   //-----------------------------------------------------------------
1529   //
1530   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1531   Double_t y,z;
1532   if (!t.GetProlongation(x,y,z)) {
1533     t.fRemoval = 10;
1534     return 0;
1535   }
1536   //
1537   //
1538   if (TMath::Abs(y)>ymax){
1539     return 0;
1540
1541     if (y > ymax) {
1542       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1543       if (!t.Rotate(fSectors->GetAlpha()))
1544         return 0;
1545     } else if (y <-ymax) {
1546       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1547       if (!t.Rotate(-fSectors->GetAlpha()))
1548         return 0;
1549     }
1550     if (!t.PropagateTo(x)) {
1551       return 0;
1552     }
1553     t.GetProlongation(x,y,z);
1554   }
1555   //
1556   // update current shape info every 3 pad-row
1557   if ( (nr%6==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1558     //    t.fCurrentSigmaY = GetSigmaY(&t);
1559     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1560     GetShape(&t,nr);
1561   }
1562   //
1563   AliTPCclusterMI *cl=0;
1564   UInt_t index=0;
1565
1566
1567   //Int_t nr2 = nr;
1568   const AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1569   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1570   Double_t  roady  =1.;
1571   Double_t  roadz = 1.;
1572   //
1573   Int_t row = nr;
1574   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1575     t.fInDead = kTRUE;
1576     t.SetClusterIndex2(row,-1);
1577     return 0;
1578   }
1579   else
1580     {
1581       if (TMath::Abs(z)>(1.05*x+10)) t.SetClusterIndex2(row,-1);
1582     }
1583   //calculate
1584
1585   if ((cl==0)&&(krow)) {
1586     //    cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1587     cl = krow.FindNearest2(y,z,roady,roadz,index);
1588
1589     if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1590   }
1591
1592   if (cl) {
1593     t.fCurrentCluster = cl;
1594     //    Int_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1595     //if (accept<3){
1596       t.SetClusterIndex2(row,index);
1597       t.fClusterPointer[row] = cl;
1598       //}
1599   }
1600   return 1;
1601 }
1602
1603
1604
1605 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
1606   //-----------------------------------------------------------------
1607   // This function tries to find a track prolongation to next pad row
1608   //-----------------------------------------------------------------
1609   t.fCurrentCluster  = 0;
1610   t.fCurrentClusterIndex1 = 0;
1611
1612   Double_t xt=t.GetX();
1613   Int_t     row = GetRowNumber(xt)-1;
1614   Double_t  ymax= GetMaxY(nr);
1615
1616   if (row < nr) return 1; // don't prolongate if not information until now -
1617   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1618     t.fRemoval =10;
1619     return 0;  // not prolongate strongly inclined tracks
1620   }
1621   if (TMath::Abs(t.GetSnp())>0.95) {
1622     t.fRemoval =10;
1623     return 0;  // not prolongate strongly inclined tracks
1624   }
1625
1626   Double_t x= GetXrow(nr);
1627   Double_t y,z;
1628   //t.PropagateTo(x+0.02);
1629   //t.PropagateTo(x+0.01);
1630   if (!t.PropagateTo(x)){
1631     return 0;
1632   }
1633   //
1634   y=t.GetY();
1635   z=t.GetZ();
1636
1637   if (TMath::Abs(y)>ymax){
1638     if (y > ymax) {
1639       t.fRelativeSector= (t.fRelativeSector+1) % fN;
1640       if (!t.Rotate(fSectors->GetAlpha()))
1641         return 0;
1642     } else if (y <-ymax) {
1643       t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
1644       if (!t.Rotate(-fSectors->GetAlpha()))
1645         return 0;
1646     }
1647     if (!t.PropagateTo(x)){
1648       return 0;
1649     }
1650     y = t.GetY();
1651   }
1652   //
1653
1654   AliTPCRow &krow=GetRow(t.fRelativeSector,nr);
1655
1656   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
1657     t.fInDead = kTRUE;
1658     t.SetClusterIndex2(nr,-1);
1659     return 0;
1660   }
1661   else
1662     {
1663       if (TMath::Abs(t.GetZ())<(1.05*t.GetX()+10)) t.fNFoundable++;
1664       else
1665         return 0;
1666     }
1667
1668   // update current
1669   if ( (nr%6==0) || t.GetNumberOfClusters()<2){
1670     //    t.fCurrentSigmaY = GetSigmaY(&t);
1671     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1672     GetShape(&t,nr);
1673   }
1674
1675   AliTPCclusterMI *cl=0;
1676   UInt_t index=0;
1677   //
1678   Double_t roady = 1.;
1679   Double_t roadz = 1.;
1680   //
1681   if (krow) {
1682     //cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1683     cl = krow.FindNearest2(y,z,roady,roadz,index);
1684   }
1685   t.fCurrentCluster  = cl;
1686   if (cl) t.fCurrentClusterIndex1 = krow.GetIndex(index);
1687   return 1;
1688 }
1689
1690
1691 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1692   //-----------------------------------------------------------------
1693   // This function tries to find a track prolongation to next pad row
1694   //-----------------------------------------------------------------
1695
1696   //update error according neighborhoud
1697
1698   if (t.fCurrentCluster) {
1699     t.fRow = nr;
1700     Bool_t accept = AcceptCluster(&t,t.fCurrentCluster,1.);
1701
1702     if (t.fCurrentCluster->IsUsed(10)){
1703       //
1704       //
1705       //  t.fErrorZ2*=2;
1706       //  t.fErrorY2*=2;
1707       t.fNShared++;
1708       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1709         t.fRemoval =10;
1710         return 0;
1711       }
1712     }
1713
1714     if (accept<3) UpdateTrack(&t,accept);
1715
1716   } else {
1717     if (fIteration==0){
1718       if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.fRemoval=10;
1719       if (  t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.fRemoval=10;
1720
1721       if (( (t.fNFoundable*0.5 > t.GetNumberOfClusters()) || t.fNoCluster>15)) t.fRemoval=10;
1722     }
1723   }
1724   return 1;
1725 }
1726
1727
1728
1729 //_____________________________________________________________________________
1730 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1731   //-----------------------------------------------------------------
1732   // This function tries to find a track prolongation.
1733   //-----------------------------------------------------------------
1734   Double_t xt=t.GetX();
1735   //
1736   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1737   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1738   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1739   //
1740   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1741
1742   Int_t first = GetRowNumber(xt)-1;
1743   for (Int_t nr= first; nr>=rf; nr-=step) {
1744     if (FollowToNext(t,nr)==0)
1745       if (!t.IsActive()) return 0;
1746
1747   }
1748   return 1;
1749 }
1750
1751
1752 //_____________________________________________________________________________
1753 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
1754   //-----------------------------------------------------------------
1755   // This function tries to find a track prolongation.
1756   //-----------------------------------------------------------------
1757   Double_t xt=t.GetX();
1758   //
1759   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1760   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1761   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1762   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1763
1764   for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
1765
1766     if (FollowToNextFast(t,nr)==0)
1767       if (!t.IsActive()) return 0;
1768
1769   }
1770   return 1;
1771 }
1772
1773
1774
1775
1776
1777 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
1778   //-----------------------------------------------------------------
1779   // This function tries to find a track prolongation.
1780   //-----------------------------------------------------------------
1781   //  Double_t xt=t.GetX();
1782   //
1783   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1784   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1785   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1786   t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN;
1787
1788   Int_t first = 0;
1789   first = t.fFirstPoint+3;
1790   //
1791   if (first<0) first=0;
1792   for (Int_t nr=first+1; nr<=rf; nr++) {
1793     //if ( (t.GetSnp()<0.9))
1794       FollowToNext(t,nr);
1795   }
1796   return 1;
1797 }
1798
1799
1800
1801
1802
1803 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
1804 {
1805   //
1806   //
1807   sum1=0;
1808   sum2=0;
1809   Int_t sum=0;
1810   //
1811   Float_t dz2 =(s1->GetZ() - s2->GetZ());
1812   dz2*=dz2;
1813
1814   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
1815   dy2*=dy2;
1816   Float_t distance = TMath::Sqrt(dz2+dy2);
1817   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
1818
1819   //  Int_t offset =0;
1820   Int_t firstpoint = TMath::Min(s1->fFirstPoint,s2->fFirstPoint);
1821   Int_t lastpoint = TMath::Max(s1->fLastPoint,s2->fLastPoint);
1822   if (lastpoint>160)
1823     lastpoint =160;
1824   if (firstpoint<0)
1825     firstpoint = 0;
1826   if (firstpoint>lastpoint) {
1827     firstpoint =lastpoint;
1828     //    lastpoint  =160;
1829   }
1830
1831
1832   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
1833     if (s1->GetClusterIndex2(i)>0) sum1++;
1834     if (s2->GetClusterIndex2(i)>0) sum2++;
1835     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
1836       sum++;
1837     }
1838   }
1839   if (sum<5) return 0;
1840
1841   Float_t summin = TMath::Min(sum1+1,sum2+1);
1842   Float_t ratio = (sum+1)/Float_t(summin);
1843   return ratio;
1844 }
1845
1846 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
1847 {
1848   //
1849   //
1850   if (TMath::Abs(s1->GetC()-s2->GetC())>0.004) return;
1851   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>0.6) return;
1852
1853   Float_t dz2 =(s1->GetZ() - s2->GetZ());
1854   dz2*=dz2;
1855   Float_t dy2 =(s1->GetY() - s2->GetY());
1856   dy2*=dy2;
1857   Float_t distance = dz2+dy2;
1858   if (distance>325.) return ; // if there are far away  - not overlap - to reduce combinatorics
1859
1860   //
1861   Int_t sumshared=0;
1862   //
1863   Int_t firstpoint = TMath::Max(s1->fFirstPoint,s2->fFirstPoint);
1864   Int_t lastpoint = TMath::Min(s1->fLastPoint,s2->fLastPoint);
1865   //
1866   if (firstpoint>=lastpoint-5) return;;
1867
1868   for (Int_t i=firstpoint;i<lastpoint;i++){
1869     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1870     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1871       sumshared++;
1872     }
1873   }
1874   if (sumshared>4){
1875     // sign clusters
1876     //
1877     for (Int_t i=firstpoint;i<lastpoint;i++){
1878       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
1879       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
1880         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
1881         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);;
1882         if (s1->IsActive()&&s2->IsActive()){
1883           p1->fIsShared = kTRUE;
1884           p2->fIsShared = kTRUE;
1885         }
1886       }
1887     }
1888   }
1889   //
1890   if (sumshared>10){
1891     for (Int_t i=0;i<4;i++){
1892       if (s1->fOverlapLabels[3*i]==0){
1893         s1->fOverlapLabels[3*i] = s2->GetLabel();
1894         s1->fOverlapLabels[3*i+1] = sumshared;
1895         s1->fOverlapLabels[3*i+2] = s2->GetUniqueID();
1896         break;
1897       }
1898     }
1899     for (Int_t i=0;i<4;i++){
1900       if (s2->fOverlapLabels[3*i]==0){
1901         s2->fOverlapLabels[3*i] = s1->GetLabel();
1902         s2->fOverlapLabels[3*i+1] = sumshared;
1903         s2->fOverlapLabels[3*i+2] = s1->GetUniqueID();
1904         break;
1905       }
1906     }
1907   }
1908
1909 }
1910
1911 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
1912 {
1913   //
1914   //sort trackss according sectors
1915   //
1916   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1917     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1918     if (!pt) continue;
1919     //if (pt) RotateToLocal(pt);
1920     pt->fSort = 0;
1921   }
1922   arr->UnSort();
1923   arr->Sort();  // sorting according z
1924   arr->Expand(arr->GetEntries());
1925   //
1926   //
1927   Int_t nseed=arr->GetEntriesFast();
1928   for (Int_t i=0; i<nseed; i++) {
1929     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1930     if (!pt) continue;
1931     for (Int_t j=0;j<=12;j++){
1932       pt->fOverlapLabels[j] =0;
1933     }
1934   }
1935   for (Int_t i=0; i<nseed; i++) {
1936     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1937     if (!pt) continue;
1938     if (pt->fRemoval>10) continue;
1939     for (Int_t j=i+1; j<nseed; j++){
1940       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1941       //      if (pt2){
1942       if (pt2->fRemoval<=10) {
1943         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
1944         SignShared(pt,pt2);
1945       }
1946     }
1947   }
1948 }
1949
1950 void  AliTPCtrackerMI::RemoveDouble(TObjArray * arr, Float_t factor1, Float_t factor2,  Int_t removalindex)
1951 {
1952   //
1953   //sort trackss according sectors
1954   //
1955   if (fDebug&1) {
1956     printf("Number of tracks before double removal- %d\n",arr->GetEntries());
1957   }
1958   //
1959   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
1960     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1961     if (!pt) continue;
1962     pt->fSort = 0;
1963   }
1964   arr->UnSort();
1965   arr->Sort();  // sorting according z
1966   arr->Expand(arr->GetEntries());
1967   //
1968   //reset overlap labels
1969   //
1970   Int_t nseed=arr->GetEntriesFast();
1971   for (Int_t i=0; i<nseed; i++) {
1972     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1973     if (!pt) continue;
1974     pt->SetUniqueID(i);
1975     for (Int_t j=0;j<=12;j++){
1976       pt->fOverlapLabels[j] =0;
1977     }
1978   }
1979   //
1980   //sign shared tracks
1981   for (Int_t i=0; i<nseed; i++) {
1982     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
1983     if (!pt) continue;
1984     if (pt->fRemoval>10) continue;
1985     Float_t deltac = pt->GetC()*0.1;
1986     for (Int_t j=i+1; j<nseed; j++){
1987       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
1988       //      if (pt2){
1989       if (pt2->fRemoval<=10) {
1990         if ( TMath::Abs(pt->fRelativeSector-pt2->fRelativeSector)>0) break;
1991         if (TMath::Abs(pt->GetC()  -pt2->GetC())>deltac) continue;
1992         if (TMath::Abs(pt->GetTgl()-pt2->GetTgl())>0.05) continue;
1993         //
1994         SignShared(pt,pt2);
1995       }
1996     }
1997   }
1998   //
1999   // remove highly shared tracks
2000   for (Int_t i=0; i<nseed; i++) {
2001     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2002     if (!pt) continue;
2003     if (pt->fRemoval>10) continue;
2004     //
2005     Int_t sumshared =0;
2006     for (Int_t j=0;j<4;j++){
2007       sumshared = pt->fOverlapLabels[j*3+1];
2008     }
2009     Float_t factor = factor1;
2010     if (pt->fRemoval>0) factor = factor2;
2011     if (sumshared/pt->GetNumberOfClusters()>factor){
2012       for (Int_t j=0;j<4;j++){
2013         if (pt->fOverlapLabels[3*j]==0) continue;
2014         if (pt->fOverlapLabels[3*j+1]<5) continue;
2015         if (pt->fRemoval==removalindex) continue;
2016         AliTPCseed * pt2 = (AliTPCseed*)arr->UncheckedAt(pt->fOverlapLabels[3*j+2]);
2017         if (!pt2) continue;
2018         if (pt2->GetSigma2C()<pt->GetSigma2C()){
2019           //      pt->fRemoval = removalindex;
2020           delete arr->RemoveAt(i);
2021           break;
2022         }
2023       }
2024     }
2025   }
2026   arr->Compress();
2027   if (fDebug&1) {
2028     printf("Number of tracks after double removal- %d\n",arr->GetEntries());
2029   }
2030 }
2031
2032
2033
2034
2035
2036
2037 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode)
2038 {
2039   //
2040   //sort tracks in array according mode criteria
2041   Int_t nseed = arr->GetEntriesFast();
2042   for (Int_t i=0; i<nseed; i++) {
2043     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2044     if (!pt) {
2045       continue;
2046     }
2047     pt->fSort = mode;
2048   }
2049   arr->UnSort();
2050   arr->Sort();
2051 }
2052
2053 void AliTPCtrackerMI::RemoveUsed(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t removalindex)
2054 {
2055
2056   //Loop over all tracks and remove "overlaps"
2057   //
2058   //
2059   Int_t nseed = arr->GetEntriesFast();
2060   Int_t good =0;
2061
2062   for (Int_t i=0; i<nseed; i++) {
2063     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2064     if (!pt) {
2065       delete arr->RemoveAt(i);
2066     }
2067     else{
2068       pt->fSort =1;
2069       pt->fBSigned = kFALSE;
2070     }
2071   }
2072   arr->Compress();
2073   nseed = arr->GetEntriesFast();
2074   arr->UnSort();
2075   arr->Sort();
2076   //
2077   //unsign used
2078   UnsignClusters();
2079   //
2080   for (Int_t i=0; i<nseed; i++) {
2081     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2082     if (!pt) {
2083       continue;
2084     }
2085     Int_t found,foundable,shared;
2086     if (pt->IsActive())
2087       pt->GetClusterStatistic(0,160,found, foundable,shared,kFALSE);
2088     else
2089       pt->GetClusterStatistic(0,160,found, foundable,shared,kTRUE);
2090     //
2091     Double_t factor = factor2;
2092     if (pt->fBConstrain) factor = factor1;
2093
2094     if ((Float_t(shared)/Float_t(found))>factor){
2095       pt->Desactivate(removalindex);
2096       continue;
2097     }
2098
2099     good++;
2100     for (Int_t i=0; i<160; i++) {
2101       Int_t index=pt->GetClusterIndex2(i);
2102       if (index<0 || index&0x8000 ) continue;
2103       AliTPCclusterMI *c= pt->fClusterPointer[i];
2104       if (!c) continue;
2105       //      if (!c->IsUsed(10)) c->Use(10);
2106       //if (pt->IsActive())
2107       c->Use(10);
2108       //else
2109       //        c->Use(5);
2110     }
2111
2112   }
2113   fNtracks = good;
2114
2115   printf("\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2116 }
2117
2118 void AliTPCtrackerMI::UnsignClusters()
2119 {
2120   //
2121   // loop over all clusters and unsign them
2122   //
2123
2124   for (Int_t sec=0;sec<fkNIS;sec++){
2125     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2126       AliTPCclusterMI *cl = fInnerSec[sec][row].fClusters1;
2127       for (Int_t icl =0;icl< fInnerSec[sec][row].fN1;icl++)
2128         //      if (cl[icl].IsUsed(10))
2129         cl[icl].Use(-1);
2130       cl = fInnerSec[sec][row].fClusters2;
2131       for (Int_t icl =0;icl< fInnerSec[sec][row].fN2;icl++)
2132         //if (cl[icl].IsUsed(10))
2133           cl[icl].Use(-1);
2134     }
2135   }
2136
2137   for (Int_t sec=0;sec<fkNOS;sec++){
2138     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2139       AliTPCclusterMI *cl = fOuterSec[sec][row].fClusters1;
2140       for (Int_t icl =0;icl< fOuterSec[sec][row].fN1;icl++)
2141         //if (cl[icl].IsUsed(10))
2142           cl[icl].Use(-1);
2143       cl = fOuterSec[sec][row].fClusters2;
2144       for (Int_t icl =0;icl< fOuterSec[sec][row].fN2;icl++)
2145         //if (cl[icl].IsUsed(10))
2146         cl[icl].Use(-1);
2147     }
2148   }
2149
2150 }
2151
2152
2153
2154 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2155 {
2156   //
2157   //sign clusters to be "used"
2158   //
2159   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2160   // loop over "primaries"
2161
2162   Float_t sumdens=0;
2163   Float_t sumdens2=0;
2164   Float_t sumn   =0;
2165   Float_t sumn2  =0;
2166   Float_t sumchi =0;
2167   Float_t sumchi2 =0;
2168
2169   Float_t sum    =0;
2170
2171   TStopwatch timer;
2172   timer.Start();
2173
2174   Int_t nseed = arr->GetEntriesFast();
2175   for (Int_t i=0; i<nseed; i++) {
2176     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2177     if (!pt) {
2178       continue;
2179     }
2180     if (!(pt->IsActive())) continue;
2181     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2182     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2183       sumdens += dens;
2184       sumdens2+= dens*dens;
2185       sumn    += pt->GetNumberOfClusters();
2186       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2187       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2188       if (chi2>5) chi2=5;
2189       sumchi  +=chi2;
2190       sumchi2 +=chi2*chi2;
2191       sum++;
2192     }
2193   }
2194
2195   Float_t mdensity = 0.9;
2196   Float_t meann    = 130;
2197   Float_t meanchi  = 1;
2198   Float_t sdensity = 0.1;
2199   Float_t smeann    = 10;
2200   Float_t smeanchi  =0.4;
2201
2202
2203   if (sum>20){
2204     mdensity = sumdens/sum;
2205     meann    = sumn/sum;
2206     meanchi  = sumchi/sum;
2207     //
2208     sdensity = sumdens2/sum-mdensity*mdensity;
2209     sdensity = TMath::Sqrt(sdensity);
2210     //
2211     smeann   = sumn2/sum-meann*meann;
2212     smeann   = TMath::Sqrt(smeann);
2213     //
2214     smeanchi = sumchi2/sum - meanchi*meanchi;
2215     smeanchi = TMath::Sqrt(smeanchi);
2216   }
2217
2218
2219   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2220   //
2221   for (Int_t i=0; i<nseed; i++) {
2222     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2223     if (!pt) {
2224       continue;
2225     }
2226     if (pt->fBSigned) continue;
2227     if (pt->fBConstrain) continue;
2228     //if (!(pt->IsActive())) continue;
2229     /*
2230     Int_t found,foundable,shared;
2231     pt->GetClusterStatistic(0,160,found, foundable,shared);
2232     if (shared/float(found)>0.3) {
2233       if (shared/float(found)>0.9 ){
2234         //delete arr->RemoveAt(i);
2235       }
2236       continue;
2237     }
2238     */
2239     Bool_t isok =kFALSE;
2240     if ( (pt->fNShared/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2241       isok = kTRUE;
2242     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->fNShared/pt->GetNumberOfClusters()<0.7))
2243       isok =kTRUE;
2244     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2245       isok =kTRUE;
2246     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2247       isok =kTRUE;
2248
2249     if (isok)
2250       for (Int_t i=0; i<160; i++) {
2251         Int_t index=pt->GetClusterIndex2(i);
2252         if (index<0) continue;
2253         AliTPCclusterMI *c= pt->fClusterPointer[i];
2254         if (!c) continue;
2255         //if (!(c->IsUsed(10))) c->Use();
2256         c->Use(10);
2257       }
2258   }
2259
2260
2261   //
2262   Double_t maxchi  = meanchi+2.*smeanchi;
2263
2264   for (Int_t i=0; i<nseed; i++) {
2265     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2266     if (!pt) {
2267       continue;
2268     }
2269     //if (!(pt->IsActive())) continue;
2270     if (pt->fBSigned) continue;
2271     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2272     if (chi>maxchi) continue;
2273
2274     Float_t bfactor=1;
2275     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->fNFoundable);
2276
2277     //sign only tracks with enoug big density at the beginning
2278
2279     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2280
2281
2282     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2283     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2284
2285     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2286     if ( (pt->fRemoval==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2287       minn=0;
2288
2289     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2290       //Int_t noc=pt->GetNumberOfClusters();
2291       pt->fBSigned = kTRUE;
2292       for (Int_t i=0; i<160; i++) {
2293
2294         Int_t index=pt->GetClusterIndex2(i);
2295         if (index<0) continue;
2296         AliTPCclusterMI *c= pt->fClusterPointer[i];
2297         if (!c) continue;
2298         //      if (!(c->IsUsed(10))) c->Use();
2299         c->Use(10);
2300       }
2301     }
2302   }
2303   //  gLastCheck = nseed;
2304   //  arr->Compress();
2305   if (fDebug>0){
2306     timer.Print();
2307   }
2308 }
2309
2310
2311 void  AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2)
2312 {
2313   // stop not active tracks
2314   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2315   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2316   Int_t nseed = arr->GetEntriesFast();
2317   //
2318   for (Int_t i=0; i<nseed; i++) {
2319     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2320     if (!pt) {
2321       continue;
2322     }
2323     if (!(pt->IsActive())) continue;
2324     StopNotActive(pt,row0,th0, th1,th2);
2325   }
2326 }
2327
2328
2329
2330 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2331  Float_t th2)
2332 {
2333   // stop not active tracks
2334   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2335   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2336   Int_t sumgood1  = 0;
2337   Int_t sumgood2  = 0;
2338   Int_t foundable = 0;
2339   Int_t maxindex = seed->fLastPoint;  //last foundable row
2340   if (seed->fNFoundable*th0 > seed->GetNumberOfClusters()) {
2341     seed->Desactivate(10) ;
2342     return;
2343   }
2344
2345   for (Int_t i=row0; i<maxindex; i++){
2346     Int_t index = seed->GetClusterIndex2(i);
2347     if (index!=-1) foundable++;
2348     //if (!c) continue;
2349     if (foundable<=30) sumgood1++;
2350     if (foundable<=50) {
2351       sumgood2++;
2352     }
2353     else{
2354       break;
2355     }
2356   }
2357   if (foundable>=30.){
2358      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2359   }
2360   if (foundable>=50)
2361     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2362 }
2363
2364
2365
2366 Int_t AliTPCtrackerMI::PropagateBack(AliESD *event)
2367 {
2368   //
2369   // back propagation of ESD tracks
2370   //
2371
2372   fEvent = event;
2373   ReadSeeds(event);
2374   PropagateBack(fSeeds);
2375   Int_t nseed = fSeeds->GetEntriesFast();
2376   for (Int_t i=0;i<nseed;i++){
2377     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2378     AliESDtrack *esd=event->GetTrack(i);
2379     seed->CookdEdx(0.02,0.06);
2380     CookLabel(seed,0.1); //For comparison only
2381     esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2382   }
2383   fEvent =0;
2384   WriteTracks();
2385   return 0;
2386 }
2387
2388
2389 void AliTPCtrackerMI::DeleteSeeds()
2390 {
2391   Int_t nseed = fSeeds->GetEntriesFast();
2392   for (Int_t i=0;i<nseed;i++){
2393     AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2394     if (seed) delete fSeeds->RemoveAt(i);
2395   }
2396   delete fSeeds;
2397   fSeeds =0;
2398 }
2399
2400 void AliTPCtrackerMI::ReadSeeds(AliESD *event)
2401 {
2402   //
2403   //read seeds from the event
2404
2405   Int_t nentr=event->GetNumberOfTracks();
2406   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
2407   if (fSeeds)
2408     DeleteSeeds();
2409   if (!fSeeds){
2410     fSeeds = new TObjArray;
2411   }
2412
2413   //  Int_t ntrk=0;
2414   for (Int_t i=0; i<nentr; i++) {
2415     AliESDtrack *esd=event->GetTrack(i);
2416     ULong_t status=esd->GetStatus();
2417     const AliTPCtrack t(*esd);
2418     AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2419     if (status==AliESDtrack::kTPCin) seed->Modify(0.8);
2420     //
2421     //
2422     // rotate to the local coordinate system
2423
2424     fSectors=fInnerSec; fN=fkNIS;
2425
2426     Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2427     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2428     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2429     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2430     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2431     alpha-=seed->GetAlpha();
2432     if (!seed->Rotate(alpha)) continue;
2433     //
2434     seed->PropagateTo(fSectors->GetX(0));
2435     //
2436     //    Int_t index = esd->GetTPCindex();
2437     //AliTPCseed * seed2= (AliTPCseed*)fSeeds->At(index);
2438
2439     fSeeds->AddLast(seed);
2440   }
2441 }
2442
2443
2444
2445 //_____________________________________________________________________________
2446 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
2447                                  Float_t deltay, Int_t ddsec) {
2448   //-----------------------------------------------------------------
2449   // This function creates track seeds.
2450   // SEEDING WITH VERTEX CONSTRAIN
2451   //-----------------------------------------------------------------
2452   // cuts[0]   - fP4 cut
2453   // cuts[1]   - tan(phi)  cut
2454   // cuts[2]   - zvertex cut
2455   // cuts[3]   - fP3 cut
2456   Int_t nin0  = 0;
2457   Int_t nin1  = 0;
2458   Int_t nin2  = 0;
2459   Int_t nin   = 0;
2460   Int_t nout1 = 0;
2461   Int_t nout2 = 0;
2462
2463   Double_t x[5], c[15];
2464   //  Int_t di = i1-i2;
2465   //
2466   AliTPCseed * seed = new AliTPCseed;
2467   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2468   Double_t cs=cos(alpha), sn=sin(alpha);
2469   //
2470   //  Double_t x1 =fOuterSec->GetX(i1);
2471   //Double_t xx2=fOuterSec->GetX(i2);
2472
2473   Double_t x1 =GetXrow(i1);
2474   Double_t xx2=GetXrow(i2);
2475
2476   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2477
2478   Int_t imiddle = (i2+i1)/2;    //middle pad row index
2479   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2480   const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
2481   //
2482   Int_t ns =sec;
2483
2484   const AliTPCRow& kr1=GetRow(ns,i1);
2485   Double_t ymax  = GetMaxY(i1)-kr1.fDeadZone-1.5;
2486   Double_t ymaxm = GetMaxY(imiddle)-kr1.fDeadZone-1.5;
2487
2488   //
2489   // change cut on curvature if it can't reach this layer
2490   // maximal curvature set to reach it
2491   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2492   if (dvertexmax*0.5*cuts[0]>0.85){
2493     cuts[0] = 0.85/(dvertexmax*0.5+1.);
2494   }
2495   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
2496
2497   //  Int_t ddsec = 1;
2498   if (deltay>0) ddsec = 0;
2499   // loop over clusters
2500   for (Int_t is=0; is < kr1; is++) {
2501     //
2502     if (kr1[is]->IsUsed(10)) continue;
2503     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2504     //if (TMath::Abs(y1)>ymax) continue;
2505
2506     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2507
2508     // find possible directions
2509     Float_t anglez = (z1-z3)/(x1-x3);
2510     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
2511     //
2512     //
2513     //find   rotation angles relative to line given by vertex and point 1
2514     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2515     Double_t dvertex  = TMath::Sqrt(dvertex2);
2516     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
2517     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);
2518
2519     //
2520     // loop over 2 sectors
2521     Int_t dsec1=-ddsec;
2522     Int_t dsec2= ddsec;
2523     if (y1<0)  dsec2= 0;
2524     if (y1>0)  dsec1= 0;
2525
2526     Double_t dddz1=0;  // direction of delta inclination in z axis
2527     Double_t dddz2=0;
2528     if ( (z1-z3)>0)
2529       dddz1 =1;
2530     else
2531       dddz2 =1;
2532     //
2533     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2534       Int_t sec2 = sec + dsec;
2535       //
2536       //      AliTPCRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2537       //AliTPCRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2538       AliTPCRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
2539       AliTPCRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2540       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2541       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2542
2543       // rotation angles to p1-p3
2544       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2545       Double_t x2,   y2,   z2;
2546       //
2547       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2548
2549       //
2550       Double_t dxx0 =  (xx2-x3)*cs13r;
2551       Double_t dyy0 =  (xx2-x3)*sn13r;
2552       for (Int_t js=index1; js < index2; js++) {
2553         const AliTPCclusterMI *kcl = kr2[js];
2554         if (kcl->IsUsed(10)) continue;
2555         //
2556         //calcutate parameters
2557         //
2558         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
2559         // stright track
2560         if (TMath::Abs(yy0)<0.000001) continue;
2561         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
2562         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2563         Double_t r02 = (0.25+y0*y0)*dvertex2;
2564         //curvature (radius) cut
2565         if (r02<r2min) continue;
2566
2567         nin0++;
2568         //
2569         Double_t c0  = 1/TMath::Sqrt(r02);
2570         if (yy0>0) c0*=-1.;
2571
2572
2573         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
2574         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2575         Double_t dfi0   = 2.*TPCFastMath::FastAsin(dvertex*c0*0.5);
2576         Double_t dfi1   = 2.*TPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2577         //
2578         //
2579         Double_t z0  =  kcl->GetZ();
2580         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
2581         if (TMath::Abs(zzzz2-z0)>0.5) continue;
2582         nin1++;
2583         //
2584         Double_t dip    = (z1-z0)*c0/dfi1;
2585         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
2586         //
2587         y2 = kcl->GetY();
2588         if (dsec==0){
2589           x2 = xx2;
2590           z2 = kcl->GetZ();
2591         }
2592         else
2593           {
2594             // rotation
2595             z2 = kcl->GetZ();
2596             x2= xx2*cs-y2*sn*dsec;
2597             y2=+xx2*sn*dsec+y2*cs;
2598           }
2599
2600         x[0] = y1;
2601         x[1] = z1;
2602         x[2] = x0;
2603         x[3] = dip;
2604         x[4] = c0;
2605         //
2606         //
2607         // do we have cluster at the middle ?
2608         Double_t ym,zm;
2609         GetProlongation(x1,xm,x,ym,zm);
2610         UInt_t dummy;
2611         AliTPCclusterMI * cm=0;
2612         if (TMath::Abs(ym)-ymaxm<0){
2613           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
2614           if ((!cm) || (cm->IsUsed(10))) {
2615             continue;
2616           }
2617         }
2618         else{
2619           // rotate y1 to system 0
2620           // get state vector in rotated system
2621           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
2622           Double_t xr2  =  x0*cs+yr1*sn*dsec;
2623           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
2624           //
2625           GetProlongation(xx2,xm,xr,ym,zm);
2626           if (TMath::Abs(ym)-ymaxm<0){
2627             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
2628             if ((!cm) || (cm->IsUsed(10))) {
2629               continue;
2630             }
2631           }
2632         }
2633
2634
2635         Double_t dym = 0;
2636         Double_t dzm = 0;
2637         if (cm){
2638           dym = ym - cm->GetY();
2639           dzm = zm - cm->GetZ();
2640         }
2641         nin2++;
2642
2643
2644         //
2645         //
2646         Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
2647         Double_t sy2=kcl->GetSigmaY2()*2.,     sz2=kcl->GetSigmaZ2()*2.;
2648         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
2649         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
2650         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
2651
2652         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2653         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2654         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2655         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2656         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2657         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2658
2659         Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2660         Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2661         Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2662         Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2663
2664         c[0]=sy1;
2665         c[1]=0.;       c[2]=sz1;
2666         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2667         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
2668                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2669         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2670         c[13]=f30*sy1*f40+f32*sy2*f42;
2671         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2672
2673         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2674
2675         UInt_t index=kr1.GetIndex(is);
2676         AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, ns*alpha+shift);
2677
2678         track->fIsSeeding = kTRUE;
2679         track->fSeed1 = i1;
2680         track->fSeed2 = i2;
2681         track->fSeedType=3;
2682
2683
2684         //if (dsec==0) {
2685           FollowProlongation(*track, (i1+i2)/2,1);
2686           Int_t foundable,found,shared;
2687           track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
2688           if ((found<0.55*foundable)  || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
2689             seed->Reset();
2690             seed->~AliTPCseed();
2691             continue;
2692           }
2693           //}
2694
2695         nin++;
2696         FollowProlongation(*track, i2,1);
2697
2698
2699         //Int_t rc = 1;
2700         track->fBConstrain =1;
2701         //      track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
2702         track->fLastPoint = i1;  // first cluster in track position
2703         track->fFirstPoint = track->fLastPoint;
2704
2705         if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2706             track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
2707             track->fNShared>0.4*track->GetNumberOfClusters() ) {
2708           seed->Reset();
2709           seed->~AliTPCseed();
2710           continue;
2711         }
2712         nout1++;
2713         // Z VERTEX CONDITION
2714         Double_t zv;
2715         zv = track->GetZ()+track->GetTgl()/track->GetC()*
2716           ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2717         if (TMath::Abs(zv-z3)>cuts[2]) {
2718           FollowProlongation(*track, TMath::Max(i2-20,0));
2719           zv = track->GetZ()+track->GetTgl()/track->GetC()*
2720             ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2721           if (TMath::Abs(zv-z3)>cuts[2]){
2722             FollowProlongation(*track, TMath::Max(i2-40,0));
2723             zv = track->GetZ()+track->GetTgl()/track->GetC()*
2724               ( asin(-track->GetEta()) - asin(track->GetX()*track->GetC()-track->GetEta()));
2725             if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->fNFoundable*0.7)){
2726               // make seed without constrain
2727               AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
2728               FollowProlongation(*track2, i2,1);
2729               track2->fBConstrain = kFALSE;
2730               track2->fSeedType = 1;
2731               arr->AddLast(track2);
2732               seed->Reset();
2733               seed->~AliTPCseed();
2734               continue;
2735             }
2736             else{
2737               seed->Reset();
2738               seed->~AliTPCseed();
2739               continue;
2740
2741             }
2742           }
2743         }
2744
2745         track->fSeedType =0;
2746         arr->AddLast(track);
2747         seed = new AliTPCseed;
2748         nout2++;
2749         // don't consider other combinations
2750         if (track->GetNumberOfClusters() > track->fNFoundable*0.8)
2751           break;
2752       }
2753     }
2754   }
2755   if (fDebug>1){
2756     //    printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
2757   }
2758   delete seed;
2759 }
2760
2761
2762 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
2763                                  Float_t deltay) {
2764
2765
2766
2767   //-----------------------------------------------------------------
2768   // This function creates track seeds.
2769   //-----------------------------------------------------------------
2770   // cuts[0]   - fP4 cut
2771   // cuts[1]   - tan(phi)  cut
2772   // cuts[2]   - zvertex cut
2773   // cuts[3]   - fP3 cut
2774
2775
2776   Int_t nin0  = 0;
2777   Int_t nin1  = 0;
2778   Int_t nin2  = 0;
2779   Int_t nin   = 0;
2780   Int_t nout1 = 0;
2781   Int_t nout2 = 0;
2782   Int_t nout3 =0;
2783   Double_t x[5], c[15];
2784   //
2785   // make temporary seed
2786   AliTPCseed * seed = new AliTPCseed;
2787   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
2788   //  Double_t cs=cos(alpha), sn=sin(alpha);
2789   //
2790   //
2791
2792   // first 3 padrows
2793   Double_t x1 = GetXrow(i1-1);
2794   const    AliTPCRow& kr1=GetRow(sec,i1-1);
2795   Double_t y1max  = GetMaxY(i1-1)-kr1.fDeadZone-1.5;
2796   //
2797   Double_t x1p = GetXrow(i1);
2798   const    AliTPCRow& kr1p=GetRow(sec,i1);
2799   //
2800   Double_t x1m = GetXrow(i1-2);
2801   const    AliTPCRow& kr1m=GetRow(sec,i1-2);
2802
2803   //
2804   //last 3 padrow for seeding
2805   AliTPCRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
2806   Double_t    x3   =  GetXrow(i1-7);
2807   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
2808   //
2809   AliTPCRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
2810   Double_t    x3p   = GetXrow(i1-6);
2811   //
2812   AliTPCRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
2813   Double_t    x3m   = GetXrow(i1-8);
2814
2815   //
2816   //
2817   // middle padrow
2818   Int_t im = i1-4;                           //middle pad row index
2819   Double_t xm         = GetXrow(im);         // radius of middle pad-row
2820   const AliTPCRow& krm=GetRow(sec,im);   //middle pad -row
2821   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
2822   //
2823   //
2824   Double_t deltax  = x1-x3;
2825   Double_t dymax   = deltax*cuts[1];
2826   Double_t dzmax   = deltax*cuts[3];
2827   //
2828   // loop over clusters
2829   for (Int_t is=0; is < kr1; is++) {
2830     //
2831     if (kr1[is]->IsUsed(10)) continue;
2832     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2833     //
2834     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2835     //
2836     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
2837     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
2838     //
2839     Double_t y3,   z3;
2840     //
2841     //
2842     UInt_t index;
2843     for (Int_t js=index1; js < index2; js++) {
2844       const AliTPCclusterMI *kcl = kr3[js];
2845       if (kcl->IsUsed(10)) continue;
2846       y3 = kcl->GetY();
2847       // apply angular cuts
2848       if (TMath::Abs(y1-y3)>dymax) continue;
2849       x3 = x3;
2850       z3 = kcl->GetZ();
2851       if (TMath::Abs(z1-z3)>dzmax) continue;
2852       //
2853       Double_t angley = (y1-y3)/(x1-x3);
2854       Double_t anglez = (z1-z3)/(x1-x3);
2855       //
2856       Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
2857       Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
2858       //
2859       Double_t yyym = angley*(xm-x1)+y1;
2860       Double_t zzzm = anglez*(xm-x1)+z1;
2861
2862       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
2863       if (!kcm) continue;
2864       if (kcm->IsUsed(10)) continue;
2865
2866       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
2867       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
2868       //
2869       //
2870       //
2871       Int_t used  =0;
2872       Int_t found =0;
2873       //
2874       // look around first
2875       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
2876                                                       anglez*(x1m-x1)+z1,
2877                                                       erry,errz,index);
2878       //
2879       if (kc1m){
2880         found++;
2881         if (kc1m->IsUsed(10)) used++;
2882       }
2883       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
2884                                                       anglez*(x1p-x1)+z1,
2885                                                       erry,errz,index);
2886       //
2887       if (kc1p){
2888         found++;
2889         if (kc1p->IsUsed(10)) used++;
2890       }
2891       if (used>1)  continue;
2892       if (found<1) continue;
2893
2894       //
2895       // look around last
2896       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
2897                                                       anglez*(x3m-x3)+z3,
2898                                                       erry,errz,index);
2899       //
2900       if (kc3m){
2901         found++;
2902         if (kc3m->IsUsed(10)) used++;
2903       }
2904       else
2905         continue;
2906       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
2907                                                       anglez*(x3p-x3)+z3,
2908                                                       erry,errz,index);
2909       //
2910       if (kc3p){
2911         found++;
2912         if (kc3p->IsUsed(10)) used++;
2913       }
2914       else
2915         continue;
2916       if (used>1)  continue;
2917       if (found<3) continue;
2918       //
2919       Double_t x2,y2,z2;
2920       x2 = xm;
2921       y2 = kcm->GetY();
2922       z2 = kcm->GetZ();
2923       //
2924
2925       x[0]=y1;
2926       x[1]=z1;
2927       x[4]=f1(x1,y1,x2,y2,x3,y3);
2928       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
2929       nin0++;
2930       //
2931       x[2]=f2(x1,y1,x2,y2,x3,y3);
2932       nin1++;
2933       //
2934       x[3]=f3n(x1,y1,x2,y2,z1,z2,x[4]);
2935       //if (TMath::Abs(x[3]) > cuts[3]) continue;
2936       nin2++;
2937       //
2938       //
2939       Double_t sy1=0.1,  sz1=0.1;
2940       Double_t sy2=0.1,  sz2=0.1;
2941       Double_t sy3=0.1,  sy=0.1, sz=0.1;
2942
2943       Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
2944       Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
2945       Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
2946       Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
2947       Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
2948       Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
2949
2950       Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
2951       Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
2952       Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
2953       Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
2954
2955       c[0]=sy1;
2956       c[1]=0.;       c[2]=sz1;
2957       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
2958       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
2959       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
2960       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
2961       c[13]=f30*sy1*f40+f32*sy2*f42;
2962       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
2963
2964       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
2965
2966       UInt_t index=kr1.GetIndex(is);
2967       AliTPCseed *track=new(seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
2968
2969       track->fIsSeeding = kTRUE;
2970
2971       nin++;
2972       FollowProlongation(*track, i1-7,1);
2973       if (track->GetNumberOfClusters() < track->fNFoundable*0.75 ||
2974           track->fNShared>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
2975         seed->Reset();
2976         seed->~AliTPCseed();
2977         continue;
2978       }
2979       nout1++;
2980       nout2++;
2981       //Int_t rc = 1;
2982       FollowProlongation(*track, i2,1);
2983       track->fBConstrain =0;
2984       track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
2985       track->fFirstPoint = track->fLastPoint;
2986
2987       if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
2988           track->GetNumberOfClusters()<track->fNFoundable*0.7 ||
2989           track->fNShared>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
2990         seed->Reset();
2991         seed->~AliTPCseed();
2992         continue;
2993       }
2994
2995       {
2996         FollowProlongation(*track, TMath::Max(i2-10,0),1);
2997         AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
2998         FollowProlongation(*track2, i2,1);
2999         track2->fBConstrain = kFALSE;
3000         track2->fSeedType = 4;
3001         arr->AddLast(track2);
3002         seed->Reset();
3003         seed->~AliTPCseed();
3004       }
3005
3006
3007       //arr->AddLast(track);
3008       //seed = new AliTPCseed;
3009       nout3++;
3010     }
3011   }
3012
3013   if (fDebug>1){
3014     //    printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3015   }
3016   delete seed;
3017 }
3018
3019
3020 //_____________________________________________________________________________
3021 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3022                                  Float_t deltay, Bool_t /*bconstrain*/) {
3023   //-----------------------------------------------------------------
3024   // This function creates track seeds - without vertex constraint
3025   //-----------------------------------------------------------------
3026   // cuts[0]   - fP4 cut        - not applied
3027   // cuts[1]   - tan(phi)  cut
3028   // cuts[2]   - zvertex cut    - not applied
3029   // cuts[3]   - fP3 cut
3030   Int_t nin0=0;
3031   Int_t nin1=0;
3032   Int_t nin2=0;
3033   Int_t nin3=0;
3034   //  Int_t nin4=0;
3035   //Int_t nin5=0;
3036
3037
3038
3039   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3040   //  Double_t cs=cos(alpha), sn=sin(alpha);
3041   Int_t row0 = (i1+i2)/2;
3042   Int_t drow = (i1-i2)/2;
3043   const AliTPCRow& kr0=fSectors[sec][row0];
3044   AliTPCRow * kr=0;
3045
3046   AliTPCpolyTrack polytrack;
3047   Int_t nclusters=fSectors[sec][row0];
3048   AliTPCseed * seed = new AliTPCseed;
3049
3050   Int_t sumused=0;
3051   Int_t cused=0;
3052   Int_t cnused=0;
3053   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3054     Int_t nfound =0;
3055     Int_t nfoundable =0;
3056     for (Int_t iter =1; iter<2; iter++){   //iterations
3057       const AliTPCRow& krm=fSectors[sec][row0-iter];
3058       const AliTPCRow& krp=fSectors[sec][row0+iter];
3059       const AliTPCclusterMI * cl= kr0[is];
3060
3061       if (cl->IsUsed(10)) {
3062         cused++;
3063       }
3064       else{
3065         cnused++;
3066       }
3067       Double_t x = kr0.GetX();
3068       // Initialization of the polytrack
3069       nfound =0;
3070       nfoundable =0;
3071       polytrack.Reset();
3072       //
3073       Double_t y0= cl->GetY();
3074       Double_t z0= cl->GetZ();
3075       Float_t erry = 0;
3076       Float_t errz = 0;
3077
3078       Double_t ymax = fSectors->GetMaxY(row0)-kr0.fDeadZone-1.5;
3079       if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;  // seed only at the edge
3080
3081       erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3082       errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3083       polytrack.AddPoint(x,y0,z0,erry, errz);
3084
3085       sumused=0;
3086       if (cl->IsUsed(10)) sumused++;
3087
3088
3089       Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3090       Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3091       //
3092       x = krm.GetX();
3093       AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3094       if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3095         erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3096         errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3097         if (cl1->IsUsed(10))  sumused++;
3098         polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3099       }
3100       //
3101       x = krp.GetX();
3102       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3103       if (cl2) {
3104         erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3105         errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3106         if (cl2->IsUsed(10)) sumused++;
3107         polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3108       }
3109       //
3110       if (sumused>0) continue;
3111       nin0++;
3112       polytrack.UpdateParameters();
3113       // follow polytrack
3114       roadz = 1.2;
3115       roady = 1.2;
3116       //
3117       Double_t yn,zn;
3118       nfoundable = polytrack.GetN();
3119       nfound     = nfoundable;
3120       //
3121       for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3122         Float_t maxdist = 0.8*(1.+3./(ddrow));
3123         for (Int_t delta = -1;delta<=1;delta+=2){
3124           Int_t row = row0+ddrow*delta;
3125           kr = &(fSectors[sec][row]);
3126           Double_t xn = kr->GetX();
3127           Double_t ymax = fSectors->GetMaxY(row)-kr->fDeadZone-1.5;
3128           polytrack.GetFitPoint(xn,yn,zn);
3129           if (TMath::Abs(yn)>ymax) continue;
3130           nfoundable++;
3131           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3132           if (cln) {
3133             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3134             if (dist<maxdist){
3135               /*
3136               erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3137               errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3138               if (cln->IsUsed(10)) {
3139                 //      printf("used\n");
3140                 sumused++;
3141                 erry*=2;
3142                 errz*=2;
3143               }
3144               */
3145               erry=0.1;
3146               errz=0.1;
3147               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3148               nfound++;
3149             }
3150           }
3151         }
3152         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;
3153         polytrack.UpdateParameters();
3154       }
3155     }
3156     if ( (sumused>3) || (sumused>0.5*nfound))  {
3157       //printf("sumused   %d\n",sumused);
3158       continue;
3159     }
3160     nin1++;
3161     Double_t dy,dz;
3162     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3163     AliTPCpolyTrack track2;
3164
3165     polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3166     if (track2.GetN()<0.5*nfoundable) continue;
3167     nin2++;
3168
3169     if ((nfound>0.6*nfoundable) &&( nfoundable>0.4*(i1-i2))) {
3170       //
3171       // test seed with and without constrain
3172       for (Int_t constrain=0; constrain<=0;constrain++){
3173         // add polytrack candidate
3174
3175         Double_t x[5], c[15];
3176         Double_t x1,x2,x3,y1,y2,y3,z1,z2,z3;
3177         track2.GetBoundaries(x3,x1);
3178         x2 = (x1+x3)/2.;
3179         track2.GetFitPoint(x1,y1,z1);
3180         track2.GetFitPoint(x2,y2,z2);
3181         track2.GetFitPoint(x3,y3,z3);
3182         //
3183         //is track pointing to the vertex ?
3184         Double_t x0,y0,z0;
3185         x0=0;
3186         polytrack.GetFitPoint(x0,y0,z0);
3187
3188         if (constrain) {
3189           x2 = x3;
3190           y2 = y3;
3191           z2 = z3;
3192
3193           x3 = 0;
3194           y3 = 0;
3195           z3 = 0;
3196         }
3197         x[0]=y1;
3198         x[1]=z1;
3199         x[4]=f1(x1,y1,x2,y2,x3,y3);
3200
3201         //      if (TMath::Abs(x[4]) >= cuts[0]) continue;  //
3202         x[2]=f2(x1,y1,x2,y2,x3,y3);
3203
3204         //if (TMath::Abs(x[4]*x1-x[2]) >= cuts[1]) continue;
3205         //x[3]=f3(x1,y1,x2,y2,z1,z2);
3206         x[3]=f3n(x1,y1,x3,y3,z1,z3,x[4]);
3207         //if (TMath::Abs(x[3]) > cuts[3]) continue;
3208
3209
3210         Double_t sy =0.1, sz =0.1;
3211         Double_t sy1=0.02, sz1=0.02;
3212         Double_t sy2=0.02, sz2=0.02;
3213         Double_t sy3=0.02;
3214
3215         if (constrain){
3216           sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3217         }
3218
3219         Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3220         Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3221         Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3222         Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3223         Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3224         Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3225
3226         Double_t f30=(f3(x1,y1+sy,x3,y3,z1,z3)-x[3])/sy;
3227         Double_t f31=(f3(x1,y1,x3,y3,z1+sz,z3)-x[3])/sz;
3228         Double_t f32=(f3(x1,y1,x3,y3+sy,z1,z3)-x[3])/sy;
3229         Double_t f34=(f3(x1,y1,x3,y3,z1,z3+sz)-x[3])/sz;
3230
3231
3232         c[0]=sy1;
3233         c[1]=0.;       c[2]=sz1;
3234         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3235         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3236         c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3237         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3238         c[13]=f30*sy1*f40+f32*sy2*f42;
3239         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3240
3241         //Int_t row1 = fSectors->GetRowNumber(x1);
3242         Int_t row1 = GetRowNumber(x1);
3243
3244         UInt_t index=0;
3245         //kr0.GetIndex(is);
3246         AliTPCseed *track=new (seed) AliTPCseed(index, x, c, x1, sec*alpha+shift);
3247         track->fIsSeeding = kTRUE;
3248         Int_t rc=FollowProlongation(*track, i2);
3249         if (constrain) track->fBConstrain =1;
3250         else
3251           track->fBConstrain =0;
3252         track->fLastPoint = row1+fInnerSec->GetNRows();  // first cluster in track position
3253         track->fFirstPoint = track->fLastPoint;
3254
3255         if (rc==0 || track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3256             track->GetNumberOfClusters() < track->fNFoundable*0.6 ||
3257             track->fNShared>0.4*track->GetNumberOfClusters()) {
3258           //delete track;
3259           seed->Reset();
3260           seed->~AliTPCseed();
3261         }
3262         else {
3263           arr->AddLast(track);
3264           seed = new AliTPCseed;
3265         }
3266         nin3++;
3267       }
3268     }  // if accepted seed
3269   }
3270   if (fDebug>1){
3271     printf("\nSeeding statiistic:\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin3);
3272   }
3273   delete seed;
3274 }
3275
3276
3277 AliTPCseed *AliTPCtrackerMI::MakeSeed(AliTPCseed *track, Float_t r0, Float_t r1, Float_t r2)
3278 {
3279   //
3280   //
3281   //reseed
3282   Int_t p0 = int(r0*track->GetNumberOfClusters());     // point 0
3283   Int_t p1 = int(r1*track->GetNumberOfClusters());
3284   Int_t p2 = int(r2*track->GetNumberOfClusters());   // last point
3285   Int_t pp2=0;
3286   Double_t  x0[3],x1[3],x2[3];
3287   x0[0]=-1;
3288   x0[0]=-1;
3289   x0[0]=-1;
3290
3291   // find track position at given ratio of the length
3292   Int_t  sec0, sec1, sec2;
3293   sec0=0;
3294   sec1=0;
3295   sec2=0;
3296   Int_t index=-1;
3297   Int_t clindex;
3298   for (Int_t i=0;i<160;i++){
3299     if (track->fClusterPointer[i]){
3300       index++;
3301       AliTPCTrackerPoint   *trpoint =track->GetTrackPoint(i);
3302       if ( (index<p0) || x0[0]<0 ){
3303         if (trpoint->GetX()>1){
3304           clindex = track->GetClusterIndex2(i);
3305           if (clindex>0){
3306             x0[0] = trpoint->GetX();
3307             x0[1] = trpoint->GetY();
3308             x0[2] = trpoint->GetZ();
3309             sec0  = ((clindex&0xff000000)>>24)%18;
3310           }
3311         }
3312       }
3313
3314       if ( (index<p1) &&(trpoint->GetX()>1)){
3315         clindex = track->GetClusterIndex2(i);
3316         if (clindex>0){
3317           x1[0] = trpoint->GetX();
3318           x1[1] = trpoint->GetY();
3319           x1[2] = trpoint->GetZ();
3320           sec1  = ((clindex&0xff000000)>>24)%18;
3321         }
3322       }
3323       if ( (index<p2) &&(trpoint->GetX()>1)){
3324         clindex = track->GetClusterIndex2(i);
3325         if (clindex>0){
3326           x2[0] = trpoint->GetX();
3327           x2[1] = trpoint->GetY();
3328           x2[2] = trpoint->GetZ();
3329           sec2  = ((clindex&0xff000000)>>24)%18;
3330           pp2 = i;
3331         }
3332       }
3333     }
3334   }
3335
3336   Double_t alpha, cs,sn, xx2,yy2;
3337   //
3338   alpha = (sec1-sec2)*fSectors->GetAlpha();
3339   cs = TMath::Cos(alpha);
3340   sn = TMath::Sin(alpha);
3341   xx2= x1[0]*cs-x1[1]*sn;
3342   yy2= x1[0]*sn+x1[1]*cs;
3343   x1[0] = xx2;
3344   x1[1] = yy2;
3345   //
3346   alpha = (sec0-sec2)*fSectors->GetAlpha();
3347   cs = TMath::Cos(alpha);
3348   sn = TMath::Sin(alpha);
3349   xx2= x0[0]*cs-x0[1]*sn;
3350   yy2= x0[0]*sn+x0[1]*cs;
3351   x0[0] = xx2;
3352   x0[1] = yy2;
3353   //
3354   //
3355   //
3356   Double_t x[5],c[15];
3357   //
3358   x[0]=x2[1];
3359   x[1]=x2[2];
3360   x[4]=f1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3361   //  if (x[4]>1) return 0;
3362   x[2]=f2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]);
3363   x[3]=f3n(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2],x[4]);
3364   //if (TMath::Abs(x[3]) > 2.2)  return 0;
3365   //if (TMath::Abs(x[2]) > 1.99) return 0;
3366   //
3367   Double_t sy =0.1,  sz =0.1;
3368   //
3369   Double_t sy1=0.02+track->GetSigmaY2(), sz1=0.02+track->GetSigmaZ2();
3370   Double_t sy2=0.01+track->GetSigmaY2(), sz2=0.01+track->GetSigmaZ2();
3371   Double_t sy3=0.01+track->GetSigmaY2();
3372   //
3373   Double_t f40=(f1(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[4])/sy;
3374   Double_t f42=(f1(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[4])/sy;
3375   Double_t f43=(f1(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[4])/sy;
3376   Double_t f20=(f2(x2[0],x2[1]+sy,x1[0],x1[1],x0[0],x0[1])-x[2])/sy;
3377   Double_t f22=(f2(x2[0],x2[1],x1[0],x1[1]+sy,x0[0],x0[1])-x[2])/sy;
3378   Double_t f23=(f2(x2[0],x2[1],x1[0],x1[1],x0[0],x0[1]+sy)-x[2])/sy;
3379   //
3380   Double_t f30=(f3(x2[0],x2[1]+sy,x0[0],x0[1],x2[2],x0[2])-x[3])/sy;
3381   Double_t f31=(f3(x2[0],x2[1],x0[0],x0[1],x2[2]+sz,x0[2])-x[3])/sz;
3382   Double_t f32=(f3(x2[0],x2[1],x0[0],x0[1]+sy,x2[2],x0[2])-x[3])/sy;
3383   Double_t f34=(f3(x2[0],x2[1],x0[0],x0[1],x2[2],x0[2]+sz)-x[3])/sz;
3384
3385
3386   c[0]=sy1;
3387   c[1]=0.;       c[2]=sz1;
3388   c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3389   c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3390   c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3391   c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3392   c[13]=f30*sy1*f40+f32*sy2*f42;
3393   c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3394
3395   //  Int_t row1 = fSectors->GetRowNumber(x2[0]);
3396   AliTPCseed *seed=new  AliTPCseed(0, x, c, x2[0], sec2*fSectors->GetAlpha()+fSectors->GetAlphaShift());
3397   //  Double_t y0,z0,y1,z1, y2,z2;
3398   //seed->GetProlongation(x0[0],y0,z0);
3399   // seed->GetProlongation(x1[0],y1,z1);
3400   //seed->GetProlongation(x2[0],y2,z2);
3401   //  seed =0;
3402   seed->fLastPoint  = pp2;
3403   seed->fFirstPoint = pp2;
3404
3405
3406   return seed;
3407 }
3408
3409 Int_t  AliTPCtrackerMI::CheckKinkPoint(AliTPCseed*seed, Float_t th)
3410 {
3411   //
3412   //
3413   //
3414   for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3415   //
3416   if (TMath::Abs(seed->GetC())>0.01) return 0;
3417   //
3418
3419   Float_t x[160], y[160], erry[160], z[160], errz[160];
3420   Int_t sec[160];
3421   Float_t xt[160], yt[160], zt[160];
3422   Int_t i1 = 200;
3423   Int_t i2 = 0;
3424   Int_t secm   = -1;
3425   Int_t padm   = -1;
3426   Int_t middle = seed->GetNumberOfClusters()/2;
3427   //
3428   //
3429   // find central sector, get local cooordinates
3430   Int_t count = 0;
3431   for (Int_t i=seed->fFirstPoint;i<=seed->fLastPoint;i++) {
3432     sec[i]= seed->GetClusterSector(i)%18;
3433     x[i]  = GetXrow(i);
3434     if (sec[i]>=0) {
3435       AliTPCclusterMI * cl = seed->fClusterPointer[i];
3436       //      if (cl==0)        cl = GetClusterMI(seed->GetClusterIndex2(i));
3437       if (cl==0) {
3438         sec[i] = -1;
3439         continue;
3440       }
3441       //
3442       //
3443       if (i>i2)  i2 = i;  //last  point with cluster
3444       if (i2<i1) i1 = i;  //first point with cluster
3445       y[i] = cl->GetY();
3446       z[i] = cl->GetZ();
3447       AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3448       xt[i] = x[i];
3449       yt[i] = point->GetY();
3450       zt[i] = point->GetZ();
3451
3452       if (point->GetX()>0){
3453         erry[i] = point->GetErrY();
3454         errz[i] = point->GetErrZ();
3455       }
3456
3457       count++;
3458       if (count<middle) {
3459         secm = sec[i];  //central sector
3460         padm = i;       //middle point with cluster
3461       }
3462     }
3463   }
3464   //
3465   // rotate position to global coordinate system connected to  sector at last the point
3466   //
3467   for (Int_t i=i1;i<=i2;i++){
3468     //
3469     if (sec[i]<0) continue;
3470     Double_t alpha = (sec[i2]-sec[i])*fSectors->GetAlpha();
3471     Double_t cs = TMath::Cos(alpha);
3472     Double_t sn = TMath::Sin(alpha);
3473     Float_t xx2= x[i]*cs+y[i]*sn;
3474     Float_t yy2= -x[i]*sn+y[i]*cs;
3475     x[i] = xx2;
3476     y[i] = yy2;
3477     //
3478     xx2= xt[i]*cs+yt[i]*sn;
3479     yy2= -xt[i]*sn+yt[i]*cs;
3480     xt[i] = xx2;
3481     yt[i] = yy2;
3482
3483   }
3484   //get "state" vector
3485   Double_t xh[5],xm = x[padm];
3486   xh[0]=yt[i2];
3487   xh[1]=zt[i2];
3488   xh[4]=f1(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3489   xh[2]=f2(xt[i2],yt[i2],xt[padm],yt[padm],xt[i1],yt[i1]);
3490   xh[3]=f3n(xt[i2],yt[i2],xt[i1],yt[i1],zt[i2],zt[i1],xh[4]);
3491   //
3492   //
3493   for (Int_t i=i1;i<=i2;i++){
3494     Double_t yy,zz;
3495     if (sec[i]<0) continue;
3496     GetProlongation(x[i2], x[i],xh,yy,zz);
3497     if (TMath::Abs(y[i]-yy)>4||TMath::Abs(z[i]-zz)>4){
3498       //Double_t xxh[5];
3499       //xxh[4]=f1old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3500       //xxh[2]=f2old(x[i2],y[i2],x[padm],y[padm],x[i1],y[i1]);
3501       printf("problem\n");
3502     }
3503     y[i] = y[i] - yy;
3504     z[i] = z[i] - zz;
3505   }
3506   Float_t dyup[160],dydown[160], dzup[160], dzdown[160];
3507   Float_t yup[160], ydown[160],  zup[160],  zdown[160];
3508
3509   AliTPCpolyTrack ptrack1,ptrack2;
3510   //
3511   // derivation up
3512   for (Int_t i=i1;i<=i2;i++){
3513     AliTPCclusterMI * cl = seed->fClusterPointer[i];
3514     if (!cl) continue;
3515     if (cl->GetType()<0) continue;
3516     if (cl->GetType()>10) continue;
3517
3518     if (sec[i]>=0){
3519       ptrack1.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3520     }
3521     if (ptrack1.GetN()>4.){
3522       ptrack1.UpdateParameters();
3523       Double_t ddy,ddz;
3524       ptrack1.GetFitDerivation(x[i]-xm,ddy,ddz);
3525       Double_t yy,zz;
3526       ptrack1.GetFitPoint(x[i]-xm,yy,zz);
3527
3528       dyup[i] = ddy;
3529       dzup[i] = ddz;
3530       yup[i]  = yy;
3531       zup[i]  = zz;
3532
3533     }
3534     else{
3535       dyup[i]=0.;  //not enough points
3536     }
3537   }
3538   //
3539   // derivation down
3540   for (Int_t i=i2;i>=i1;i--){
3541     AliTPCclusterMI * cl = seed->fClusterPointer[i];
3542     if (!cl) continue;
3543     if (cl->GetType()<0) continue;
3544     if (cl->GetType()>10) continue;
3545     if (sec[i]>=0){
3546       ptrack2.AddPoint(x[i]-xm,y[i],z[i],0.1,0.1);
3547     }
3548     if (ptrack2.GetN()>4){
3549       ptrack2.UpdateParameters();
3550       Double_t ddy,ddz;
3551       ptrack2.GetFitDerivation(x[i]-xm,ddy,ddz);
3552       Double_t yy,zz;
3553       ptrack2.GetFitPoint(x[i]-xm,yy,zz);
3554
3555       dydown[i] = ddy;
3556       dzdown[i] = ddz;
3557       ydown[i]  = yy;
3558       zdown[i]  = zz;
3559     }
3560     else{
3561       dydown[i]=0.;  //not enough points
3562     }
3563   }
3564   //
3565   //
3566   // find maximal difference of the derivation
3567   for (Int_t i=0;i<12;i++) seed->fKinkPoint[i]=0;
3568
3569
3570   for (Int_t i=i1+10;i<i2-10;i++){
3571     if ( (TMath::Abs(dydown[i])<0.00000001)  ||  (TMath::Abs(dyup[i])<0.00000001) ||i<30)continue;
3572     //    printf("%f\t%f\t%f\t%f\t%f\n",x[i],dydown[i],dyup[i],dzdown[i],dzup[i]);
3573     //
3574     Float_t ddy = TMath::Abs(dydown[i]-dyup[i]);
3575     Float_t ddz = TMath::Abs(dzdown[i]-dzup[i]);
3576     if ( (ddy+ddz)> th){
3577       seed->fKinkPoint[0] = i;
3578       seed->fKinkPoint[1] = ddy;
3579       seed->fKinkPoint[2] = ddz;
3580       th = ddy+ddz;
3581     }
3582   }
3583
3584   if (fTreeDebug){
3585     //
3586     //write information to the debug tree
3587     TBranch * br = fTreeDebug->GetBranch("debug");
3588     TClonesArray * arr = new TClonesArray("AliTPCTrackPoint2");
3589     arr->ExpandCreateFast(i2-i1);
3590     br->SetAddress(&arr);
3591     //
3592     AliTPCclusterMI cldummy;
3593     cldummy.SetQ(0);
3594     AliTPCTrackPoint2 pdummy;
3595     pdummy.GetTPoint().fIsShared = 10;
3596     //
3597     Double_t alpha = sec[i2]*fSectors->GetAlpha();
3598     Double_t cs    = TMath::Cos(alpha);
3599     Double_t sn    = TMath::Sin(alpha);
3600
3601     for (Int_t i=i1;i<i2;i++){
3602       AliTPCTrackPoint2 *trpoint = (AliTPCTrackPoint2*)arr->UncheckedAt(i-i1);
3603       //cluster info
3604       AliTPCclusterMI * cl0 = seed->fClusterPointer[i];
3605       //
3606       AliTPCTrackerPoint * point = seed->GetTrackPoint(i);
3607
3608       if (cl0){
3609         Double_t x = GetXrow(i);
3610         trpoint->GetTPoint() = *point;
3611         trpoint->GetCPoint() = *cl0;
3612         trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
3613         trpoint->fID    = seed->GetUniqueID();
3614         trpoint->fLab   = seed->GetLabel();
3615         //
3616         trpoint->fGX =  cs *x + sn*point->GetY();
3617         trpoint->fGY = -sn *x + cs*point->GetY() ;
3618         trpoint->fGZ = point->GetZ();
3619         //
3620         trpoint->fDY = y[i];
3621         trpoint->fDZ = z[i];
3622         //
3623         trpoint->fDYU = dyup[i];
3624         trpoint->fDZU = dzup[i];
3625         //
3626         trpoint->fDYD = dydown[i];
3627         trpoint->fDZD = dzdown[i];
3628         //
3629         if (TMath::Abs(dyup[i])>0.00000000001 &&TMath::Abs(dydown[i])>0.00000000001){
3630           trpoint->fDDY = dydown[i]-dyup[i];
3631           trpoint->fDDZ = dzdown[i]-dzup[i];
3632         }else{
3633           trpoint->fDDY = 0.;
3634           trpoint->fDDZ = 0.;
3635         }
3636       }
3637       else{
3638         *trpoint = pdummy;
3639         trpoint->GetCPoint()= cldummy;
3640         trpoint->fID = -1;
3641       }
3642       //
3643     }
3644     fTreeDebug->Fill();
3645   }
3646
3647
3648   return 0;
3649
3650 }
3651
3652
3653
3654
3655
3656 AliTPCseed*  AliTPCtrackerMI::ReSeed(AliTPCseed *t)
3657 {
3658   //
3659   // reseed - refit -  track
3660   //
3661   Int_t first = 0;
3662   //  Int_t last  = fSectors->GetNRows()-1;
3663   //
3664   if (fSectors == fOuterSec){
3665     first = TMath::Max(first, t->fFirstPoint-fInnerSec->GetNRows());
3666     //last  =
3667   }
3668   else
3669     first = t->fFirstPoint;
3670   //
3671   AliTPCseed * seed = MakeSeed(t,0.1,0.5,0.9);
3672   FollowBackProlongation(*t,fSectors->GetNRows()-1);
3673   t->Reset(kFALSE);
3674   FollowProlongation(*t,first);
3675   return seed;
3676 }
3677
3678
3679
3680
3681
3682
3683
3684 //_____________________________________________________________________________
3685 Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
3686   //-----------------------------------------------------------------
3687   // This function reades track seeds.
3688   //-----------------------------------------------------------------
3689   TDirectory *savedir=gDirectory;
3690
3691   TFile *in=(TFile*)inp;
3692   if (!in->IsOpen()) {
3693      cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
3694      return 1;
3695   }
3696
3697   in->cd();
3698   TTree *seedTree=(TTree*)in->Get("Seeds");
3699   if (!seedTree) {
3700      cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
3701      cerr<<"can't get a tree with track seeds !\n";
3702      return 2;
3703   }
3704   AliTPCtrack *seed=new AliTPCtrack;
3705   seedTree->SetBranchAddress("tracks",&seed);
3706
3707   if (fSeeds==0) fSeeds=new TObjArray(15000);
3708
3709   Int_t n=(Int_t)seedTree->GetEntries();
3710   for (Int_t i=0; i<n; i++) {
3711      seedTree->GetEvent(i);
3712      fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
3713   }
3714
3715   delete seed;
3716   delete seedTree;
3717   savedir->cd();
3718   return 0;
3719 }
3720
3721 //_____________________________________________________________________________
3722 Int_t AliTPCtrackerMI::Clusters2Tracks() {
3723   //-----------------------------------------------------------------
3724   // This is a track finder.
3725   //-----------------------------------------------------------------
3726   TDirectory *savedir=gDirectory;
3727   TStopwatch timer;
3728   //
3729   if (!fInput) SetIO();  //set default IO using loaders
3730   if (!fInput){
3731      cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
3732      return 1;
3733   }
3734   LoadClusters();
3735   //
3736   fIteration = 0;
3737   fSeeds = Tracking();
3738
3739
3740   printf("Time for tracking: \t");timer.Print();timer.Start();
3741
3742   //activate again some tracks
3743   for (Int_t i=0; i<fSeeds->GetEntriesFast(); i++) {
3744     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3745     if (!pt) continue;
3746     Int_t nc=t.GetNumberOfClusters();
3747     if (nc<20) {
3748       delete fSeeds->RemoveAt(i);
3749       continue;
3750     }
3751     if (pt->fRemoval==10) {
3752       if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
3753         pt->Desactivate(10);  // make track again active
3754       else{
3755         pt->Desactivate(20);
3756         delete fSeeds->RemoveAt(i);
3757       }
3758     }
3759   }
3760   RemoveDouble(fSeeds,0.2,0.6,11);
3761   //RemoveUsed(fSeeds,0.9,0.9,6);
3762   //RemoveUsed(fSeeds,0.8,0.8,6);
3763   //RemoveUsed(fSeeds,0.7,0.7,6);
3764   RemoveUsed(fSeeds,0.5,0.5,6);
3765
3766   //
3767   Int_t nseed=fSeeds->GetEntriesFast();
3768   Int_t found = 0;
3769   for (Int_t i=0; i<nseed; i++) {
3770     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3771     if (!pt) continue;
3772     Int_t nc=t.GetNumberOfClusters();
3773     if (nc<15) {
3774       delete fSeeds->RemoveAt(i);
3775       continue;
3776     }
3777     CookLabel(pt,0.1); //For comparison only
3778     //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3779     if ((pt->IsActive() || (pt->fRemoval==10) )){
3780       cerr<<found++<<'\r';
3781     }
3782     else
3783       delete fSeeds->RemoveAt(i);
3784     pt->fLab2 = i;
3785   }
3786
3787
3788   //RemoveOverlap(fSeeds,0.99,7,kTRUE);
3789   SignShared(fSeeds);
3790   //RemoveUsed(fSeeds,0.9,0.9,6);
3791   //
3792   nseed=fSeeds->GetEntriesFast();
3793   found = 0;
3794   for (Int_t i=0; i<nseed; i++) {
3795     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3796     if (!pt) continue;
3797     Int_t nc=t.GetNumberOfClusters();
3798     if (nc<15) {
3799       delete fSeeds->RemoveAt(i);
3800       continue;
3801     }
3802     t.SetUniqueID(i);
3803     t.CookdEdx(0.02,0.6);
3804     //    CheckKinkPoint(&t,0.05);
3805     //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3806     if ((pt->IsActive() || (pt->fRemoval==10) )){
3807       cerr<<found++<<'\r';
3808     }
3809     else
3810       delete fSeeds->RemoveAt(i);
3811     pt->fLab2 = i;
3812   }
3813
3814   SortTracks(fSeeds, 1);
3815
3816   /*
3817   fIteration = 1;
3818   PrepareForBackProlongation(fSeeds,0.5);
3819   PropagateBack(fSeeds);
3820   printf("Time for back propagation: \t");timer.Print();timer.Start();
3821
3822   fIteration = 2;
3823
3824   PrepareForProlongation(fSeeds,1.);
3825   PropagateForward();
3826
3827   fSectors = fOuterSec;
3828   ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
3829   fSectors = fInnerSec;
3830   ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
3831   printf("Time for FORWARD propagation: \t");timer.Print();timer.Start();
3832   // RemoveUsed(fSeeds,0.7,0.7,6);
3833   //RemoveOverlap(fSeeds,0.9,7,kTRUE);
3834
3835   nseed=fSeeds->GetEntriesFast();
3836   found = 0;
3837   for (Int_t i=0; i<nseed; i++) {
3838     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;
3839     if (!pt) continue;
3840     Int_t nc=t.GetNumberOfClusters();
3841     if (nc<15) {
3842       delete fSeeds->RemoveAt(i);
3843       continue;
3844     }
3845     t.CookdEdx(0.02,0.6);
3846     //    CookLabel(pt,0.1); //For comparison only
3847     //if ((pt->IsActive() || (pt->fRemoval==10) )&& nc>50 &&pt->GetNumberOfClusters()>0.4*pt->fNFoundable){
3848     if ((pt->IsActive() || (pt->fRemoval==10) )){
3849       cerr<<found++<<'\r';
3850     }
3851     else
3852       delete fSeeds->RemoveAt(i);
3853     pt->fLab2 = i;
3854   }
3855   */
3856
3857   //  fNTracks = found;
3858   printf("Time for overlap removal, track writing and dedx cooking: \t"); timer.Print();timer.Start();
3859   //
3860   if (fOutput) {
3861     WriteTracks();
3862   }
3863   if (!fNewIO)  fOutput->Write();
3864   else
3865     AliRunLoader::GetDetectorLoader("TPC",AliConfig::fgkDefaultEventFolderName)->WriteTracks("OVERWRITE");
3866
3867
3868   cerr<<"Number of found tracks : "<<"\t"<<found<<endl;
3869   savedir->cd();
3870   //if (seedtree) delete seedtree;
3871   //  UnloadClusters();
3872   //printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
3873
3874   return 0;
3875 }
3876
3877 void AliTPCtrackerMI::Tracking(TObjArray * arr)
3878 {
3879   //
3880   // tracking of the seeds
3881   //
3882
3883   fSectors = fOuterSec;
3884   ParallelTracking(arr,150,63);
3885   fSectors = fOuterSec;
3886   ParallelTracking(arr,63,0);
3887 }
3888
3889 TObjArray * AliTPCtrackerMI::Tracking(Int_t seedtype, Int_t i1, Int_t i2, Float_t cuts[4], Float_t dy, Int_t dsec)
3890 {
3891   //
3892   //
3893   //tracking routine
3894   TObjArray * arr = new TObjArray;
3895   //
3896   fSectors = fOuterSec;
3897   TStopwatch timer;
3898   timer.Start();
3899   for (Int_t sec=0;sec<fkNOS;sec++){
3900     if (seedtype==3) MakeSeeds3(arr,sec,i1,i2,cuts,dy, dsec);
3901     if (seedtype==4) MakeSeeds5(arr,sec,i1,i2,cuts,dy);
3902     if (seedtype==2) MakeSeeds2(arr,sec,i1,i2,cuts,dy);
3903   }
3904   if (fDebug>0){
3905     printf("\nSeeding - %d\t%d\t%d\t%d\n",seedtype,i1,i2,arr->GetEntriesFast());
3906     timer.Print();
3907     timer.Start();
3908   }
3909   Tracking(arr);
3910   if (fDebug>0){
3911     timer.Print();
3912   }
3913
3914   return arr;
3915 }
3916
3917 TObjArray * AliTPCtrackerMI::Tracking()
3918 {
3919   //
3920   //
3921   TStopwatch timer;
3922   timer.Start();
3923   Int_t nup=fOuterSec->GetNRows()+fInnerSec->GetNRows();
3924
3925   TObjArray * seeds = new TObjArray;
3926   TObjArray * arr=0;
3927
3928   Int_t gap =20;
3929   Float_t cuts[4];
3930   cuts[0] = 0.002;
3931   cuts[1] = 1.5;
3932   cuts[2] = 3.;
3933   cuts[3] = 3.;
3934   Float_t fnumber  = 3.0;
3935   Float_t fdensity = 3.0;
3936
3937   //
3938   //find primaries
3939   cuts[0]=0.0066;
3940   for (Int_t delta = 0; delta<18; delta+=6){
3941     //
3942     cuts[0]=0.0070;
3943     cuts[1] = 1.5;
3944     arr = Tracking(3,nup-1-delta,nup-1-delta-gap,cuts,-1,1);
3945     SumTracks(seeds,arr);
3946     SignClusters(seeds,fnumber,fdensity);
3947     //
3948     for (Int_t i=2;i<6;i+=2){
3949       // seed high pt tracks
3950       cuts[0]=0.0022;
3951       cuts[1]=0.3;
3952       arr = Tracking(3,nup-i-delta,nup-i-delta-gap,cuts,-1,0);
3953       SumTracks(seeds,arr);
3954       SignClusters(seeds,fnumber,fdensity);
3955     }
3956   }
3957   fnumber  = 4;
3958   fdensity = 4.;
3959   //  RemoveUsed(seeds,0.9,0.9,1);
3960   //  UnsignClusters();
3961   //  SignClusters(seeds,fnumber,fdensity);
3962
3963   //find primaries
3964   cuts[0]=0.0077;
3965   for (Int_t delta = 20; delta<120; delta+=10){
3966     //
3967     // seed high pt tracks
3968     cuts[0]=0.0060;
3969     cuts[1]=0.3;
3970     cuts[2]=6.;
3971     arr = Tracking(3,nup-delta,nup-delta-gap,cuts,-1);
3972     SumTracks(seeds,arr);
3973     SignClusters(seeds,fnumber,fdensity);
3974
3975     cuts[0]=0.003;
3976     cuts[1]=0.3;
3977     cuts[2]=6.;
3978     arr = Tracking(3,nup-delta-5,nup-delta-5-gap,cuts,-1);
3979     SumTracks(seeds,arr);
3980     SignClusters(seeds,fnumber,fdensity);
3981   }
3982
3983   cuts[0] = 0.01;
3984   cuts[1] = 2.0;
3985   cuts[2] = 3.;
3986   cuts[3] = 2.0;
3987   fnumber  = 2.;
3988   fdensity = 2.;
3989
3990   if (fDebug>0){
3991     printf("\n\nPrimary seeding\t%d\n\n",seeds->GetEntriesFast());
3992     timer.Print();
3993     timer.Start();
3994   }
3995   //  RemoveUsed(seeds,0.75,0.75,1);
3996   //UnsignClusters();
3997   //SignClusters(seeds,fnumber,fdensity);
3998
3999   // find secondaries
4000
4001   cuts[0] = 0.3;
4002   cuts[1] = 1.5;
4003   cuts[2] = 3.;
4004   cuts[3] = 1.5;
4005
4006   arr = Tracking(4,nup-1,nup-1-gap,cuts,-1);
4007   SumTracks(seeds,arr);
4008   SignClusters(seeds,fnumber,fdensity);
4009   //
4010   arr = Tracking(4,nup-2,nup-2-gap,cuts,-1);
4011   SumTracks(seeds,arr);
4012   SignClusters(seeds,fnumber,fdensity);
4013   //
4014   arr = Tracking(4,nup-3,nup-3-gap,cuts,-1);
4015   SumTracks(seeds,arr);
4016   SignClusters(seeds,fnumber,fdensity);
4017   //
4018
4019
4020   for (Int_t delta = 3; delta<30; delta+=5){
4021     //
4022     cuts[0] = 0.3;
4023     cuts[1] = 1.5;
4024     cuts[2] = 3.;
4025     cuts[3] = 1.5;
4026     arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4027     SumTracks(seeds,arr);
4028     SignClusters(seeds,fnumber,fdensity);
4029     //
4030     arr = Tracking(4,nup-3-delta,nup-5-delta-gap,cuts,4);
4031     SumTracks(seeds,arr);
4032     SignClusters(seeds,fnumber,fdensity);
4033     //
4034   }
4035   fnumber  = 1;
4036   fdensity = 1;
4037   //
4038   // change cuts
4039   fnumber  = 2.;
4040   fdensity = 2.;
4041   cuts[0]=0.0080;
4042
4043   // find secondaries
4044   for (Int_t delta = 30; delta<70; delta+=10){
4045     //
4046     cuts[0] = 0.3;
4047     cuts[1] = 1.5;
4048     cuts[2] = 3.;
4049     cuts[3] = 1.5;
4050     arr = Tracking(4,nup-1-delta,nup-1-delta-gap,cuts,-1);
4051     SumTracks(seeds,arr);
4052     SignClusters(seeds,fnumber,fdensity);
4053     //
4054     arr = Tracking(4,nup-5-delta,nup-5-delta-gap,cuts,5 );
4055     SumTracks(seeds,arr);
4056     SignClusters(seeds,fnumber,fdensity);
4057   }
4058
4059   if (fDebug>0){
4060     printf("\n\nSecondary seeding\t%d\n\n",seeds->GetEntriesFast());
4061     timer.Print();
4062     timer.Start();
4063   }
4064
4065   return seeds;
4066   //
4067
4068 }
4069
4070
4071 void AliTPCtrackerMI::SumTracks(TObjArray *arr1,TObjArray *arr2)
4072 {
4073   //
4074   //sum tracks to common container
4075   //remove suspicious tracks
4076   Int_t nseed = arr2->GetEntriesFast();
4077   for (Int_t i=0;i<nseed;i++){
4078     AliTPCseed *pt=(AliTPCseed*)arr2->UncheckedAt(i);
4079     if (pt){
4080
4081       // NORMAL ACTIVE TRACK
4082       if (pt->IsActive()){
4083         arr1->AddLast(arr2->RemoveAt(i));
4084         continue;
4085       }
4086       //remove not usable tracks
4087       if (pt->fRemoval!=10){
4088         delete arr2->RemoveAt(i);
4089         continue;
4090       }
4091       // REMOVE VERY SHORT  TRACKS
4092       if (pt->GetNumberOfClusters()<20){
4093         delete arr2->RemoveAt(i);
4094         continue;
4095       }
4096       // ENABLE ONLY ENOUGH GOOD STOPPED TRACKS
4097       if (pt->GetDensityFirst(20)>0.8 || pt->GetDensityFirst(30)>0.8 || pt->GetDensityFirst(40)>0.7)
4098         arr1->AddLast(arr2->RemoveAt(i));
4099       else{
4100         delete arr2->RemoveAt(i);
4101       }
4102     }
4103   }
4104   delete arr2;
4105 }
4106
4107
4108
4109 void  AliTPCtrackerMI::ParallelTracking(TObjArray * arr, Int_t rfirst, Int_t rlast)
4110 {
4111   //
4112   // try to track in parralel
4113
4114   Int_t nseed=arr->GetEntriesFast();
4115   //prepare seeds for tracking
4116   for (Int_t i=0; i<nseed; i++) {
4117     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
4118     if (!pt) continue;
4119     if (!t.IsActive()) continue;
4120     // follow prolongation to the first layer
4121     if ( (fSectors ==fInnerSec) || (t.fFirstPoint-fParam->GetNRowLow()>rfirst+1) )
4122       FollowProlongation(t, rfirst+1);
4123   }
4124
4125
4126   //
4127   for (Int_t nr=rfirst; nr>=rlast; nr--){
4128     // make indexes with the cluster tracks for given
4129
4130     // find nearest cluster
4131     for (Int_t i=0; i<nseed; i++) {
4132       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt;
4133       if (!pt) continue;
4134       if (!pt->IsActive()) continue;
4135       //      if ( (fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
4136       if (pt->fRelativeSector>17) {
4137         continue;
4138       }
4139       UpdateClusters(t,nr);
4140     }
4141     // prolonagate to the nearest cluster - if founded
4142     for (Int_t i=0; i<nseed; i++) {
4143       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
4144       if (!pt) continue;
4145       if (!pt->IsActive()) continue;
4146       // if ((fSectors ==fOuterSec) && (pt->fFirstPoint-fParam->GetNRowLow())<nr) continue;
4147       if (pt->fRelativeSector>17) {
4148         continue;
4149       }
4150       FollowToNextCluster(*pt,nr);
4151     }
4152   }
4153 }
4154
4155 void AliTPCtrackerMI::PrepareForBackProlongation(TObjArray * arr,Float_t fac)
4156 {
4157   //
4158   //
4159   // if we use TPC track itself we have to "update" covariance
4160   //
4161   Int_t nseed= arr->GetEntriesFast();
4162   for (Int_t i=0;i<nseed;i++){
4163     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4164     if (pt) {
4165       pt->Modify(fac);
4166       //
4167       //rotate to current local system at first accepted  point
4168       Int_t index  = pt->GetClusterIndex2(pt->fFirstPoint);
4169       Int_t sec    = (index&0xff000000)>>24;
4170       sec = sec%18;
4171       Float_t angle1 = fInnerSec->GetAlpha()*sec+fInnerSec->GetAlphaShift();
4172       if (angle1>TMath::Pi())
4173         angle1-=2.*TMath::Pi();
4174       Float_t angle2 = pt->GetAlpha();
4175
4176       if (TMath::Abs(angle1-angle2)>0.001){
4177         pt->Rotate(angle1-angle2);
4178         //angle2 = pt->GetAlpha();
4179         //pt->fRelativeSector = pt->GetAlpha()/fInnerSec->GetAlpha();
4180         //if (pt->GetAlpha()<0)
4181         //  pt->fRelativeSector+=18;
4182         //sec = pt->fRelativeSector;
4183       }
4184
4185     }
4186
4187   }
4188
4189
4190 }
4191 void AliTPCtrackerMI::PrepareForProlongation(TObjArray * arr, Float_t fac)
4192 {
4193   //
4194   //
4195   // if we use TPC track itself we have to "update" covariance
4196   //
4197   Int_t nseed= arr->GetEntriesFast();
4198   for (Int_t i=0;i<nseed;i++){
4199     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4200     if (pt) {
4201       pt->Modify(fac);
4202       pt->fFirstPoint = pt->fLastPoint;
4203     }
4204
4205   }
4206
4207
4208 }
4209
4210 Int_t AliTPCtrackerMI::PropagateBack(TObjArray * arr)
4211 {
4212   //
4213   // make back propagation
4214   //
4215   Int_t nseed= arr->GetEntriesFast();
4216   for (Int_t i=0;i<nseed;i++){
4217     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4218     if (pt) {
4219       AliTPCseed *pt2 = new AliTPCseed(*pt);
4220       fSectors = fInnerSec;
4221       FollowBackProlongation(*pt,fSectors->GetNRows()-1);
4222       fSectors = fOuterSec;
4223       FollowBackProlongation(*pt,fSectors->GetNRows()-1);
4224       fSectors = fOuterSec;
4225       if (pt->GetNumberOfClusters()<35 && pt->GetLabel()>0 ){
4226         printf("\n%d",pt->GetLabel());
4227         fSectors = fInnerSec;
4228         FollowBackProlongation(*pt2,fSectors->GetNRows()-1);
4229         fSectors = fOuterSec;
4230         FollowBackProlongation(*pt2,fSectors->GetNRows()-1);
4231         fSectors = fOuterSec;
4232       }
4233     }
4234   }
4235   return 0;
4236 }
4237
4238
4239 Int_t AliTPCtrackerMI::PropagateForward2(TObjArray * arr)
4240 {
4241   //
4242   // make forward propagation
4243   //
4244   Int_t nseed= arr->GetEntriesFast();
4245   for (Int_t i=0;i<nseed;i++){
4246     AliTPCseed *pt = (AliTPCseed*)arr->UncheckedAt(i);
4247     if (pt) {
4248       AliTPCseed *pt2 = new AliTPCseed(*pt);
4249       fSectors = fOuterSec;
4250       FollowProlongation(*pt,0);
4251       fSectors = fOuterSec;
4252       FollowProlongation(*pt,0);
4253       fSectors = fInnerSec;
4254       if (pt->GetNumberOfClusters()<35 && pt->GetLabel()>0 ){
4255         printf("\n%d",pt->GetLabel());
4256         fSectors = fOuterSec;
4257         FollowProlongation(*pt2,0);
4258         fSectors = fOuterSec;
4259         FollowProlongation(*pt2,0);
4260         fSectors = fOuterSec;
4261       }
4262     }
4263   }
4264   return 0;
4265 }
4266
4267
4268 Int_t AliTPCtrackerMI::PropagateForward()
4269 {
4270   fSectors = fOuterSec;
4271   ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
4272   fSectors = fInnerSec;
4273   ParallelTracking(fSeeds,fSectors->GetNRows()-1,0);
4274   //WriteTracks();
4275   return 1;
4276 }
4277
4278
4279
4280
4281
4282
4283 Int_t AliTPCtrackerMI::PropagateBack(AliTPCseed * pt, Int_t row0, Int_t row1)
4284 {
4285   //
4286   // make back propagation, in between row0 and row1
4287   //
4288
4289   if (pt) {
4290     fSectors = fInnerSec;
4291     Int_t  r1;
4292     //
4293     if (row1<fSectors->GetNRows())
4294       r1 = row1;
4295     else
4296       r1 = fSectors->GetNRows()-1;
4297
4298     if (row0<fSectors->GetNRows()&& r1>0 )
4299       FollowBackProlongation(*pt,r1);
4300     if (row1<=fSectors->GetNRows())
4301       return 0;
4302     //
4303     r1 = row1 - fSectors->GetNRows();
4304     if (r1<=0) return 0;
4305     if (r1>=fOuterSec->GetNRows()) return 0;
4306     fSectors = fOuterSec;
4307     return FollowBackProlongation(*pt,r1);
4308   }
4309   return 0;
4310 }
4311
4312
4313
4314
4315 void  AliTPCtrackerMI::GetShape(AliTPCseed * seed, Int_t row)
4316 {
4317   //
4318   //
4319   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4320   //  Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
4321   Float_t padlength =  GetPadPitchLength(row);
4322   //
4323   Float_t sresy = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4324   Float_t angulary  = seed->GetSnp();
4325   angulary = angulary*angulary/(1-angulary*angulary);
4326   seed->fCurrentSigmaY2 = sd2+padlength*padlength*angulary/12.+sresy*sresy;
4327   //
4328   Float_t sresz = fParam->GetZSigma();
4329   Float_t angularz  = seed->GetTgl();
4330   seed->fCurrentSigmaZ2 = sd2+padlength*padlength*angularz*angularz*(1+angulary)/12.+sresz*sresz;
4331   /*
4332   Float_t wy = GetSigmaY(seed);
4333   Float_t wz = GetSigmaZ(seed);
4334   wy*=wy;
4335   wz*=wz;
4336   if (TMath::Abs(wy/seed->fCurrentSigmaY2-1)>0.0001 || TMath::Abs(wz/seed->fCurrentSigmaZ2-1)>0.0001 ){
4337     printf("problem\n");
4338   }
4339   */
4340 }
4341
4342
4343 Float_t  AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
4344 {
4345   //
4346   //
4347   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4348   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
4349   Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
4350   Float_t angular  = seed->GetSnp();
4351   angular = angular*angular/(1-angular*angular);
4352   //  angular*=angular;
4353   //angular  = TMath::Sqrt(angular/(1-angular));
4354   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
4355   return res;
4356 }
4357 Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
4358 {
4359   //
4360   //
4361   Float_t sd2 = TMath::Abs((fParam->GetZLength()-TMath::Abs(seed->GetZ())))*fParam->GetDiffL()*fParam->GetDiffL();
4362   Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
4363   Float_t sres = fParam->GetZSigma();
4364   Float_t angular  = seed->GetTgl();
4365   Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
4366   return res;
4367 }
4368
4369
4370
4371
4372 //__________________________________________________________________________
4373 void AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong) const {
4374   //--------------------------------------------------------------------
4375   //This function "cooks" a track label. If label<0, this track is fake.
4376   //--------------------------------------------------------------------
4377   Int_t noc=t->GetNumberOfClusters();
4378   if (noc<10){
4379     printf("\nnot founded prolongation\n\n\n");
4380     t->Dump();
4381     return ;
4382   }
4383   Int_t lb[160];
4384   Int_t mx[160];
4385   AliTPCclusterMI *clusters[160];
4386   //
4387   for (Int_t i=0;i<160;i++) {
4388     clusters[i]=0;
4389     lb[i]=mx[i]=0;
4390   }
4391
4392   Int_t i;
4393   Int_t current=0;
4394   for (i=0; i<160 && current<noc; i++) {
4395
4396      Int_t index=t->GetClusterIndex2(i);
4397      if (index<=0) continue;
4398      if (index&0x8000) continue;
4399      //
4400      //clusters[current]=GetClusterMI(index);
4401      if (t->fClusterPointer[i]){
4402        clusters[current]=t->fClusterPointer[i];
4403        current++;
4404      }
4405   }
4406   noc = current;
4407
4408   Int_t lab=123456789;
4409   for (i=0; i<noc; i++) {
4410     AliTPCclusterMI *c=clusters[i];
4411     if (!c) continue;
4412     lab=TMath::Abs(c->GetLabel(0));
4413     Int_t j;
4414     for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
4415     lb[j]=lab;
4416     (mx[j])++;
4417   }
4418
4419   Int_t max=0;
4420   for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
4421
4422   for (i=0; i<noc; i++) {
4423     AliTPCclusterMI *c=clusters[i];
4424     if (!c) continue;
4425     if (TMath::Abs(c->GetLabel(1)) == lab ||
4426         TMath::Abs(c->GetLabel(2)) == lab ) max++;
4427   }
4428
4429   if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
4430
4431   else {
4432      Int_t tail=Int_t(0.10*noc);
4433      max=0;
4434      Int_t ind=0;
4435      for (i=1; i<=160&&ind<tail; i++) {
4436        //       AliTPCclusterMI *c=clusters[noc-i];
4437        AliTPCclusterMI *c=clusters[i];
4438        if (!c) continue;
4439        if (lab == TMath::Abs(c->GetLabel(0)) ||
4440            lab == TMath::Abs(c->GetLabel(1)) ||
4441            lab == TMath::Abs(c->GetLabel(2))) max++;
4442        ind++;
4443      }
4444      if (max < Int_t(0.5*tail)) lab=-lab;
4445   }
4446
4447   t->SetLabel(lab);
4448
4449   //  delete[] lb;
4450   //delete[] mx;
4451   //delete[] clusters;
4452 }
4453
4454 //_________________________________________________________________________
4455 void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
4456   //-----------------------------------------------------------------------
4457   // Setup inner sector
4458   //-----------------------------------------------------------------------
4459   if (f==0) {
4460      fAlpha=par->GetInnerAngle();
4461      fAlphaShift=par->GetInnerAngleShift();
4462      fPadPitchWidth=par->GetInnerPadPitchWidth();
4463      fPadPitchLength=par->GetInnerPadPitchLength();
4464      fN=par->GetNRowLow();
4465      fRow=new AliTPCRow[fN];
4466      for (Int_t i=0; i<fN; i++) {
4467        fRow[i].SetX(par->GetPadRowRadiiLow(i));
4468        fRow[i].fDeadZone =1.5;  //1.5 cm of dead zone
4469      }
4470   } else {
4471      fAlpha=par->GetOuterAngle();
4472      fAlphaShift=par->GetOuterAngleShift();
4473      fPadPitchWidth  = par->GetOuterPadPitchWidth();
4474      fPadPitchLength = par->GetOuter1PadPitchLength();
4475      f1PadPitchLength = par->GetOuter1PadPitchLength();
4476      f2PadPitchLength = par->GetOuter2PadPitchLength();
4477
4478      fN=par->GetNRowUp();
4479      fRow=new AliTPCRow[fN];
4480      for (Int_t i=0; i<fN; i++) {
4481        fRow[i].SetX(par->GetPadRowRadiiUp(i));
4482        fRow[i].fDeadZone =1.5;  // 1.5 cm of dead zone
4483      }
4484   }
4485 }
4486
4487
4488 AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
4489   //
4490 }
4491
4492
4493
4494 //_________________________________________________________________________
4495 void
4496 AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
4497   //-----------------------------------------------------------------------
4498   // Insert a cluster into this pad row in accordence with its y-coordinate
4499   //-----------------------------------------------------------------------
4500   if (fN==kMaxClusterPerRow) {
4501     cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
4502   }
4503   if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
4504   Int_t i=Find(c->GetZ());
4505   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
4506   memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
4507   fIndex[i]=index; fClusters[i]=c; fN++;
4508 }
4509
4510
4511 //___________________________________________________________________
4512 Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
4513   //-----------------------------------------------------------------------
4514   // Return the index of the nearest cluster
4515   //-----------------------------------------------------------------------
4516   if (fN==0) return 0;
4517   if (z <= fClusters[0]->GetZ()) return 0;
4518   if (z > fClusters[fN-1]->GetZ()) return fN;
4519   Int_t b=0, e=fN-1, m=(b+e)/2;
4520   for (; b<e; m=(b+e)/2) {
4521     if (z > fClusters[m]->GetZ()) b=m+1;
4522     else e=m;
4523   }
4524   return m;
4525 }
4526
4527
4528
4529 //___________________________________________________________________
4530 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
4531   //-----------------------------------------------------------------------
4532   // Return the index of the nearest cluster in z y
4533   //-----------------------------------------------------------------------
4534   Float_t maxdistance = roady*roady + roadz*roadz;
4535
4536   AliTPCclusterMI *cl =0;
4537   for (Int_t i=Find(z-roadz); i<fN; i++) {
4538       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4539       if (c->GetZ() > z+roadz) break;
4540       if ( (c->GetY()-y) >  roady ) continue;
4541       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4542       if (maxdistance>distance) {
4543         maxdistance = distance;
4544         cl=c;
4545       }
4546   }
4547   return cl;
4548 }
4549
4550 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4551 {
4552   //-----------------------------------------------------------------------
4553   // Return the index of the nearest cluster in z y
4554   //-----------------------------------------------------------------------
4555   Float_t maxdistance = roady*roady + roadz*roadz;
4556   Int_t iz1 = TMath::Max(fFastCluster[Int_t(z-roadz+254.5)]-1,0);
4557   Int_t iz2 = TMath::Min(fFastCluster[Int_t(z+roadz+255.5)]+1,fN);
4558
4559   AliTPCclusterMI *cl =0;
4560   //FindNearest3(y,z,roady,roadz,index);
4561   //  for (Int_t i=Find(z-roadz); i<fN; i++) {
4562   for (Int_t i=iz1; i<iz2; i++) {
4563       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4564       if (c->GetZ() > z+roadz) break;
4565       if ( c->GetY()-y >  roady ) continue;
4566       if ( y-c->GetY() >  roady ) continue;
4567       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4568       if (maxdistance>distance) {
4569         maxdistance = distance;
4570         cl=c;
4571         index =i;
4572         //roady = TMath::Sqrt(maxdistance);
4573       }
4574   }
4575   return cl;
4576 }
4577
4578
4579
4580 AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const
4581 {
4582   //-----------------------------------------------------------------------
4583   // Return the index of the nearest cluster in z y
4584   //-----------------------------------------------------------------------
4585   Float_t maxdistance = roady*roady + roadz*roadz;
4586   //  Int_t iz = Int_t(z+255.);
4587   AliTPCclusterMI *cl =0;
4588   for (Int_t i=Find(z-roadz); i<fN; i++) {
4589     //for (Int_t i=fFastCluster[iz-2]; i<fFastCluster[iz+2]; i++) {
4590       AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
4591       if (c->GetZ() > z+roadz) break;
4592       if ( c->GetY()-y >  roady ) continue;
4593       if ( y-c->GetY() >  roady ) continue;
4594       Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
4595       if (maxdistance>distance) {
4596         maxdistance = distance;
4597         cl=c;
4598         index =i;
4599         //roady = TMath::Sqrt(maxdistance);
4600       }
4601   }
4602   return cl;
4603 }
4604
4605
4606
4607
4608 AliTPCseed::AliTPCseed():AliTPCtrack(){
4609   //
4610   fRow=0;
4611   fRemoval =0;
4612   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
4613   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
4614
4615   fPoints = 0;
4616   fEPoints = 0;
4617   fNFoundable =0;
4618   fNShared  =0;
4619   //  fTrackPoints =0;
4620   fRemoval = 0;
4621   fSort =0;
4622   fFirstPoint =0;
4623   fNoCluster =0;
4624   fBSigned = kFALSE;
4625   fSeed1 =-1;
4626   fSeed2 =-1;
4627 }
4628
4629 AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
4630   fPoints = 0;
4631   fEPoints = 0;
4632   fNShared  =0;
4633   //  fTrackPoints =0;
4634   fRemoval =0;
4635   fSort =0;
4636   for (Int_t i=0;i<160;i++) {
4637     fClusterPointer[i] = 0;
4638     Int_t index = t.GetClusterIndex(i);
4639     if (index>0) {
4640       SetClusterIndex2(i,index);
4641     }
4642     else{
4643       SetClusterIndex2(i,-3);
4644     }
4645   }
4646   fFirstPoint =0;
4647   fNoCluster =0;
4648   fBSigned = kFALSE;
4649   fSeed1 =-1;
4650   fSeed2 =-1;
4651 }
4652
4653 AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
4654   fRow=0;
4655   for (Int_t i=0;i<160;i++) {
4656     fClusterPointer[i] = 0;
4657     Int_t index = t.GetClusterIndex(i);
4658     SetClusterIndex2(i,index);
4659   }
4660
4661   fPoints = 0;
4662   fEPoints = 0;
4663   fNFoundable =0;
4664   fNShared  =0;
4665   //  fTrackPoints =0;
4666   fRemoval =0;
4667   fSort = 0;
4668   fFirstPoint =0;
4669   fNoCluster =0;
4670   fBSigned = kFALSE;
4671   fSeed1 =-1;
4672   fSeed2 =-1;
4673 }
4674
4675 AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15],
4676                                         Double_t xr, Double_t alpha):
4677   AliTPCtrack(index, xx, cc, xr, alpha) {
4678   //
4679   //
4680   fRow =0;
4681   for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
4682   for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
4683   fPoints = 0;
4684   fEPoints = 0;
4685   fNFoundable =0;
4686   fNShared  = 0;
4687   //  fTrackPoints =0;
4688   fRemoval =0;
4689   fSort =0;
4690   fFirstPoint =0;
4691   //  fHelixIn = new TClonesArray("AliHelix",0);
4692   //fHelixOut = new TClonesArray("AliHelix",0);
4693   fNoCluster =0;
4694   fBSigned = kFALSE;
4695   fSeed1 =-1;
4696   fSeed2 =-1;
4697 }
4698
4699 AliTPCseed::~AliTPCseed(){
4700   if (fPoints) delete fPoints;
4701   fPoints =0;
4702   if (fEPoints) delete fEPoints;
4703   fEPoints = 0;
4704   fNoCluster =0;
4705 }
4706
4707 AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
4708 {
4709   //
4710   //
4711   return &fTrackPoints[i];
4712 }
4713
4714 void AliTPCseed::RebuildSeed()
4715 {
4716   //
4717   // rebuild seed to be ready for storing
4718   AliTPCclusterMI cldummy;
4719   cldummy.SetQ(0);
4720   AliTPCTrackPoint pdummy;
4721   pdummy.GetTPoint().fIsShared = 10;
4722   for (Int_t i=0;i<160;i++){
4723     AliTPCclusterMI * cl0 = fClusterPointer[i];
4724     AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
4725     if (cl0){
4726       trpoint->GetTPoint() = *(GetTrackPoint(i));
4727       trpoint->GetCPoint() = *cl0;
4728       trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
4729     }
4730     else{
4731       *trpoint = pdummy;
4732       trpoint->GetCPoint()= cldummy;
4733     }
4734
4735   }
4736
4737 }
4738
4739
4740 Double_t AliTPCseed::GetDensityFirst(Int_t n)
4741 {
4742   //
4743   //
4744   // return cluster for n rows bellow first point
4745   Int_t nfoundable = 1;
4746   Int_t nfound      = 1;
4747   for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
4748     Int_t index = GetClusterIndex2(i);
4749     if (index!=-1) nfoundable++;
4750     if (index>0) nfound++;
4751   }
4752   if (nfoundable<n) return 0;
4753   return Double_t(nfound)/Double_t(nfoundable);
4754
4755 }
4756
4757
4758 void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
4759 {
4760   // get cluster stat.  on given region
4761   //
4762   found       = 0;
4763   foundable   = 0;
4764   shared      =0;
4765   for (Int_t i=first;i<last; i++){
4766     Int_t index = GetClusterIndex2(i);
4767     if (index!=-1) foundable++;
4768     if (fClusterPointer[i]) {
4769       found++;
4770     }
4771     else
4772       continue;
4773
4774     if (fClusterPointer[i]->IsUsed(10)) {
4775       shared++;
4776       continue;
4777     }
4778     if (!plus2) continue; //take also neighborhoud
4779     //
4780     if ( (i>0) && fClusterPointer[i-1]){
4781       if (fClusterPointer[i-1]->IsUsed(10)) {
4782         shared++;
4783         continue;
4784       }
4785     }
4786     if ( fClusterPointer[i+1]){
4787       if (fClusterPointer[i+1]->IsUsed(10)) {
4788         shared++;
4789         continue;
4790       }
4791     }
4792
4793   }
4794   if (shared>found){
4795     printf("problem\n");
4796   }
4797 }
4798
4799 //_____________________________________________________________________________
4800 void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
4801   //-----------------------------------------------------------------
4802   // This funtion calculates dE/dX within the "low" and "up" cuts.
4803   //-----------------------------------------------------------------
4804
4805   Float_t amp[200];
4806   Float_t angular[200];
4807   Float_t weight[200];
4808   Int_t index[200];
4809   //Int_t nc = 0;
4810   //  TClonesArray & arr = *fPoints;
4811   Float_t meanlog = 100.;
4812
4813   Float_t mean[4]  = {0,0,0,0};
4814   Float_t sigma[4] = {1000,1000,1000,1000};
4815   Int_t nc[4]      = {0,0,0,0};
4816   Float_t norm[4]    = {1000,1000,1000,1000};
4817   //
4818   //
4819   fNShared =0;
4820
4821   for (Int_t of =0; of<4; of++){
4822     for (Int_t i=of+i1;i<i2;i+=4)
4823       {
4824         Int_t index = fIndex[i];
4825         if (index<0||index&0x8000) continue;
4826
4827         //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
4828         AliTPCTrackerPoint * point = GetTrackPoint(i);
4829         //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
4830         //AliTPCTrackerPoint * pointp = 0;
4831         //if (i<159) pointp = GetTrackPoint(i+1);
4832
4833         if (point==0) continue;
4834         AliTPCclusterMI * cl = fClusterPointer[i];
4835         if (cl==0) continue;
4836         if (onlyused && (!cl->IsUsed(10))) continue;
4837         if (cl->IsUsed(11)) {
4838           fNShared++;
4839           continue;
4840         }
4841         Int_t   type   = cl->GetType();
4842         //if (point->fIsShared){
4843         //  fNShared++;
4844         //  continue;
4845         //}
4846         //if (pointm)
4847         //  if (pointm->fIsShared) continue;
4848         //if (pointp)
4849         //  if (pointp->fIsShared) continue;
4850
4851         if (type<0) continue;
4852         //if (type>10) continue;
4853         //if (point->GetErrY()==0) continue;
4854         //if (point->GetErrZ()==0) continue;
4855
4856         //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
4857         //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
4858         //if ((ddy*ddy+ddz*ddz)>10) continue;
4859
4860
4861         //      if (point->GetCPoint().GetMax()<5) continue;
4862         if (cl->GetMax()<5) continue;
4863         Float_t angley = point->GetAngleY();
4864         Float_t anglez = point->GetAngleZ();
4865
4866         Float_t rsigmay2 =  point->GetSigmaY();
4867         Float_t rsigmaz2 =  point->GetSigmaZ();
4868         /*
4869         Float_t ns = 1.;
4870         if (pointm){
4871           rsigmay +=  pointm->GetTPoint().GetSigmaY();
4872           rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
4873           ns+=1.;
4874         }
4875         if (pointp){
4876           rsigmay +=  pointp->GetTPoint().GetSigmaY();
4877           rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
4878           ns+=1.;
4879         }
4880         rsigmay/=ns;
4881         rsigmaz/=ns;
4882         */
4883
4884         Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
4885
4886         Float_t ampc   = 0;     // normalization to the number of electrons
4887         if (i>64){
4888           //      ampc = 1.*point->GetCPoint().GetMax();
4889           ampc = 1.*cl->GetMax();
4890           //ampc = 1.*point->GetCPoint().GetQ();
4891           //      AliTPCClusterPoint & p = point->GetCPoint();
4892           //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
4893           // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
4894           //Float_t dz =
4895           //  TMath::Abs( Int_t(iz) - iz + 0.5);
4896           //ampc *= 1.15*(1-0.3*dy);
4897           //ampc *= 1.15*(1-0.3*dz);
4898           //      Float_t zfactor = (1.05-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
4899           //ampc               *=zfactor;
4900         }
4901         else{
4902           //ampc = 1.0*point->GetCPoint().GetMax();
4903           ampc = 1.0*cl->GetMax();
4904           //ampc = 1.0*point->GetCPoint().GetQ();
4905           //AliTPCClusterPoint & p = point->GetCPoint();
4906           // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
4907           //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
4908           //Float_t dz =
4909           //  TMath::Abs( Int_t(iz) - iz + 0.5);
4910
4911           //ampc *= 1.15*(1-0.3*dy);
4912           //ampc *= 1.15*(1-0.3*dz);
4913           //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
4914           //ampc               *=zfactor;
4915
4916         }
4917         ampc *= 2.0;     // put mean value to channel 50
4918         //ampc *= 0.58;     // put mean value to channel 50
4919         Float_t w      =  1.;
4920         //      if (type>0)  w =  1./(type/2.-0.5);
4921         //      Float_t z = TMath::Abs(cl->GetZ());
4922         if (i<64) {
4923           ampc /= 0.6;
4924           //ampc /= (1+0.0008*z);
4925         } else
4926           if (i>128){
4927             ampc /=1.5;
4928             //ampc /= (1+0.0008*z);
4929           }else{
4930             //ampc /= (1+0.0008*z);
4931           }
4932
4933         if (type<0) {  //amp at the border - lower weight
4934           // w*= 2.;
4935
4936           continue;
4937         }
4938         if (rsigma>1.5) ampc/=1.3;  // if big backround
4939         amp[nc[of]]        = ampc;
4940         angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
4941         weight[nc[of]]     = w;
4942         nc[of]++;
4943       }
4944
4945     TMath::Sort(nc[of],amp,index,kFALSE);
4946     Float_t sumamp=0;
4947     Float_t sumamp2=0;
4948     Float_t sumw=0;
4949     //meanlog = amp[index[Int_t(nc[of]*0.33)]];
4950     meanlog = 50;
4951     for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
4952       Float_t ampl      = amp[index[i]]/angular[index[i]];
4953       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
4954       //
4955       sumw    += weight[index[i]];
4956       sumamp  += weight[index[i]]*ampl;
4957       sumamp2 += weight[index[i]]*ampl*ampl;
4958       norm[of]    += angular[index[i]]*weight[index[i]];
4959     }
4960     if (sumw<1){
4961       SetdEdx(0);
4962     }
4963     else {
4964       norm[of] /= sumw;
4965       mean[of]  = sumamp/sumw;
4966       sigma[of] = sumamp2/sumw-mean[of]*mean[of];
4967       if (sigma[of]>0.1)
4968         sigma[of] = TMath::Sqrt(sigma[of]);
4969       else
4970         sigma[of] = 1000;
4971
4972     mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
4973     //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
4974     //mean *=(1-0.1*(norm-1.));
4975     }
4976   }
4977
4978   Float_t dedx =0;
4979   fSdEdx =0;
4980   fMAngular =0;
4981   //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
4982   //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
4983
4984
4985   //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/
4986   //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
4987
4988   Int_t norm2 = 0;
4989   Int_t norm3 = 0;
4990   for (Int_t i =0;i<4;i++){
4991     if (nc[i]>2&&nc[i]<1000){
4992       dedx      += mean[i] *nc[i];
4993       fSdEdx    += sigma[i]*(nc[i]-2);
4994       fMAngular += norm[i] *nc[i];
4995       norm2     += nc[i];
4996       norm3     += nc[i]-2;
4997     }
4998     fDEDX[i]  = mean[i];
4999     fSDEDX[i] = sigma[i];
5000     fNCDEDX[i]= nc[i];
5001   }
5002
5003   if (norm3>0){
5004     dedx   /=norm2;
5005     fSdEdx /=norm3;
5006     fMAngular/=norm2;
5007   }
5008   else{
5009     SetdEdx(0);
5010     return;
5011   }
5012   //  Float_t dedx1 =dedx;
5013   /*
5014   dedx =0;
5015   for (Int_t i =0;i<4;i++){
5016     if (nc[i]>2&&nc[i]<1000){
5017       mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
5018       dedx      += mean[i] *nc[i];
5019     }
5020     fDEDX[i]  = mean[i];
5021   }
5022   dedx /= norm2;
5023   */
5024
5025
5026   SetdEdx(dedx);
5027
5028   //mi deDX
5029
5030
5031
5032   //Very rough PID
5033   Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
5034
5035   if (p<0.6) {
5036     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5037     if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
5038     SetMass(0.93827); return;
5039   }
5040
5041   if (p<1.2) {
5042     if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
5043     SetMass(0.93827); return;
5044   }
5045
5046   SetMass(0.13957); return;
5047
5048 }
5049
5050
5051
5052 /*
5053
5054
5055
5056 void AliTPCseed::CookdEdx2(Double_t low, Double_t up) {
5057   //-----------------------------------------------------------------
5058   // This funtion calculates dE/dX within the "low" and "up" cuts.
5059   //-----------------------------------------------------------------
5060
5061   Float_t amp[200];
5062   Float_t angular[200];
5063   Float_t weight[200];
5064   Int_t index[200];
5065   Bool_t inlimit[200];
5066   for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
5067   for (Int_t i=0;i<200;i++) amp[i]=10000;
5068   for (Int_t i=0;i<200;i++) angular[i]= 1;;
5069
5070
5071   //
5072   Float_t meanlog = 100.;
5073   Int_t indexde[4]={0,64,128,160};
5074
5075   Float_t amean     =0;
5076   Float_t asigma    =0;
5077   Float_t anc       =0;
5078   Float_t anorm     =0;
5079
5080   Float_t mean[4]  = {0,0,0,0};
5081   Float_t sigma[4] = {1000,1000,1000,1000};
5082   Int_t nc[4]      = {0,0,0,0};
5083   Float_t norm[4]    = {1000,1000,1000,1000};
5084   //
5085   //
5086   fNShared =0;
5087
5088   //  for (Int_t of =0; of<3; of++){
5089   //  for (Int_t i=indexde[of];i<indexde[of+1];i++)
5090   for (Int_t i =0; i<160;i++)
5091     {
5092         AliTPCTrackPoint * point = GetTrackPoint(i);
5093         if (point==0) continue;
5094         if (point->fIsShared){
5095           fNShared++;
5096           continue;
5097         }
5098         Int_t   type   = point->GetCPoint().GetType();
5099         if (type<0) continue;
5100         if (point->GetCPoint().GetMax()<5) continue;
5101         Float_t angley = point->GetTPoint().GetAngleY();
5102         Float_t anglez = point->GetTPoint().GetAngleZ();
5103         Float_t rsigmay =  point->GetCPoint().GetSigmaY();
5104         Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
5105         Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
5106
5107         Float_t ampc   = 0;     // normalization to the number of electrons
5108         if (i>64){
5109           ampc =  point->GetCPoint().GetMax();
5110         }
5111         else{
5112           ampc = point->GetCPoint().GetMax();
5113         }
5114         ampc *= 2.0;     // put mean value to channel 50
5115         //      ampc *= 0.565;     // put mean value to channel 50
5116
5117         Float_t w      =  1.;
5118         Float_t z = TMath::Abs(point->GetCPoint().GetZ());
5119         if (i<64) {
5120           ampc /= 0.63;
5121         } else
5122           if (i>128){
5123             ampc /=1.51;
5124           }
5125         if (type<0) {  //amp at the border - lower weight
5126           continue;
5127         }
5128         if (rsigma>1.5) ampc/=1.3;  // if big backround
5129         angular[i]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
5130         amp[i]        = ampc/angular[i];
5131         weight[i]     = w;
5132         anc++;
5133     }
5134
5135   TMath::Sort(159,amp,index,kFALSE);
5136   for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){
5137     inlimit[index[i]] = kTRUE;  // take all clusters
5138   }
5139
5140   //  meanlog = amp[index[Int_t(anc*0.3)]];
5141   meanlog =10000.;
5142   for (Int_t of =0; of<3; of++){
5143     Float_t sumamp=0;
5144     Float_t sumamp2=0;
5145     Float_t sumw=0;
5146    for (Int_t i=indexde[of];i<indexde[of+1];i++)
5147       {
5148         if (inlimit[i]==kFALSE) continue;
5149         Float_t ampl      = amp[i];
5150         ///angular[i];
5151         ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
5152         //
5153         sumw    += weight[i];
5154         sumamp  += weight[i]*ampl;
5155         sumamp2 += weight[i]*ampl*ampl;
5156         norm[of]    += angular[i]*weight[i];
5157         nc[of]++;
5158       }
5159    if (sumw<1){
5160      SetdEdx(0);
5161    }
5162    else {
5163      norm[of] /= sumw;
5164      mean[of]  = sumamp/sumw;
5165      sigma[of] = sumamp2/sumw-mean[of]*mean[of];
5166      if (sigma[of]>0.1)
5167        sigma[of] = TMath::Sqrt(sigma[of]);
5168      else
5169        sigma[of] = 1000;
5170      mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
5171    }
5172   }
5173
5174   Float_t dedx =0;
5175   fSdEdx =0;
5176   fMAngular =0;
5177   //
5178   Int_t norm2 = 0;
5179   Int_t norm3 = 0;
5180   Float_t www[3] = {12.,14.,17.};
5181   //Float_t www[3] = {1.,1.,1.};
5182
5183   for (Int_t i =0;i<3;i++){
5184     if (nc[i]>2&&nc[i]<1000){
5185       dedx      += mean[i] *nc[i]*www[i]/sigma[i];
5186       fSdEdx    += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
5187       fMAngular += norm[i] *nc[i];
5188       norm2     += nc[i]*www[i]/sigma[i];
5189       norm3     += (nc[i]-2)*www[i]/sigma[i];
5190     }
5191     fDEDX[i]  = mean[i];
5192     fSDEDX[i] = sigma[i];
5193     fNCDEDX[i]= nc[i];
5194   }
5195
5196   if (norm3>0){
5197     dedx   /=norm2;
5198     fSdEdx /=norm3;
5199     fMAngular/=norm2;
5200   }
5201   else{
5202     SetdEdx(0);
5203     return;
5204   }
5205   //  Float_t dedx1 =dedx;
5206
5207   dedx =0;
5208   Float_t norm4 = 0;
5209   for (Int_t i =0;i<3;i++){
5210     if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
5211       //mean[i]   = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
5212       dedx      += mean[i] *(nc[i])/(sigma[i]);
5213       norm4     += (nc[i])/(sigma[i]);
5214     }
5215     fDEDX[i]  = mean[i];
5216   }
5217   if (norm4>0) dedx /= norm4;
5218
5219
5220
5221   SetdEdx(dedx);
5222
5223   //mi deDX
5224
5225 }
5226
5227 */