TPC/AliTPCtrackerMI.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16
  17 //-------------------------------------------------------
  18 //          Implementation of the TPC tracker
  19 //
  20 //   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
  21 //
  22 //  AliTPC parallel tracker
  23 //
  24 //  The track fitting is based on Kalaman filtering approach
  25
  26 //  The track finding steps:
  27 //      1. Seeding - with and without vertex constraint
  28 //                 - seeding with vertex constain done at first n^2 proble
  29 //                 - seeding without vertex constraint n^3 problem
  30 //      2. Tracking - follow prolongation road - find cluster - update kalman track
  31
  32 //  The seeding and tracking is repeated several times, in different seeding region.
  33 //  This approach enables to find the track which cannot be seeded in some region of TPC
  34 //  This can happen because of low momenta (track do not reach outer radius), or track is currently in the ded region between sectors, or the track is for the moment overlapped with other track (seed quality is poor) ...
  35
  36 //  With this approach we reach almost 100 % efficiency also for high occupancy events.
  37 //  (If the seeding efficiency in a region is about 90 % than with logical or of several
  38 //  regions we will reach 100% (in theory - supposing independence)
  39
  40 //  Repeating several seeding - tracking procedures some of the tracks can be find
  41 //  several times.
  42
  43 //  The procedures to remove multi find tacks are impremented:
  44 //  RemoveUsed2 - fast procedure n problem -
  45 //                 Algorithm - Sorting tracks according quality
  46 //                             remove tracks with some shared fraction
  47 //                             Sharing in respect to all tacks
  48 //                             Signing clusters in gold region
  49 //  FindSplitted - slower algorithm n^2
  50 //                 Sort the tracks according quality
  51 //                 Loop over pair of tracks
  52 //                 If overlap with other track bigger than threshold - remove track
  53 //
  54 //  FindCurling  - Finds the pair of tracks which are curling
  55 //               - About 10% of tracks can be find with this procedure
  56 //                 The combinatorial background is too big to be used in High
  57 //                  multiplicity environment
  58 //               - n^2 problem - Slow procedure - currently it is disabled because of
  59 //                  low efficiency
  60 //
  61 //  The number of splitted tracks can be reduced disabling the sharing of the cluster.
  62 //  tpcRecoParam-> SetClusterSharing(kFALSE);
  63 //  IT IS HIGHLY non recomended to use it in high flux enviroonment
  64 //  Even using this switch some tracks can be found more than once
  65 //  (because of multiple seeding and low quality tracks which will not cross full chamber)
  66 //
  67 //
  68 // The tracker itself can be debugged  - the information about tracks can be stored in several // phases of the reconstruction
  69 // To enable storage of the TPC tracks in the ESD friend track
  70 // use AliTPCReconstructor::SetStreamLevel(n); where nis bigger 0
  71 //
  72 // The debug level -  different procedure produce tree for numerical debugging
  73 //                    To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
  74 //
  75
  76 //
  77 // Adding systematic errors to the covariance:
  78 //
  79 // The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
  80 // of the tracks (not to the clusters as they are dependent):
  81 // The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
  82 // The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
  83 // The default values are 0.
  84 //
  85 // The sytematic errors are added to the covariance matrix in following places:
  86 //
  87 // 1. During fisrt itteration - AliTPCtrackerMI::FillESD
  88 // 2. Second iteration -
  89 //      2.a ITS->TPC   - AliTPCtrackerMI::ReadSeeds
  90 //      2.b TPC->TRD   - AliTPCtrackerMI::PropagateBack
  91 // 3. Third iteration  -
  92 //      3.a TRD->TPC   - AliTPCtrackerMI::ReadSeeds
  93 //      3.b TPC->ITS   - AliTPCtrackerMI::RefitInward
  94 //
  95 // There are several places in the code which can be numerically debuged
  96 // This code is keeped in order to enable code development and to check the calibration implementtion
  97 //
  98 //      1. ErrParam stream (Log level 9) - dump information about
  99 //         1.a) cluster
 100 //         2.a) cluster error estimate
 101 //         3.a) cluster shape estimate
 102 //
 103 //
 104 //-------------------------------------------------------
 105
 106
 107 /* $Id$ */
 108
 109 #include "Riostream.h"
 110 #include <TClonesArray.h>
 111 #include <TFile.h>
 112 #include <TObjArray.h>
 113 #include <TTree.h>
 114 #include "AliLog.h"
 115 #include "AliComplexCluster.h"
 116 #include "AliESDEvent.h"
 117 #include "AliESDtrack.h"
 118 #include "AliESDVertex.h"
 119 #include "AliKink.h"
 120 #include "AliV0.h"
 121 #include "AliHelix.h"
 122 #include "AliRunLoader.h"
 123 #include "AliTPCClustersRow.h"
 124 #include "AliTPCParam.h"
 125 #include "AliTPCReconstructor.h"
 126 #include "AliTPCpolyTrack.h"
 127 #include "AliTPCreco.h"
 128 #include "AliTPCseed.h"
 129
 130 #include "AliTPCtrackerSector.h"
 131 #include "AliTPCtrackerMI.h"
 132 #include "TStopwatch.h"
 133 #include "AliTPCReconstructor.h"
 134 #include "AliPID.h"
 135 #include "TTreeStream.h"
 136 #include "AliAlignObj.h"
 137 #include "AliTrackPointArray.h"
 138 #include "TRandom.h"
 139 #include "AliTPCcalibDB.h"
 140 #include "AliTPCTransform.h"
 141 #include "AliTPCClusterParam.h"
 142
 143 //
 144
 145 ClassImp(AliTPCtrackerMI)
 146
 147
 148
 149 class AliTPCFastMath {
 150 public:
 151   AliTPCFastMath();
 152   static Double_t FastAsin(Double_t x);
 153  private:
 154   static Double_t fgFastAsin[20000];  //lookup table for fast asin computation
 155 };
 156
 157 Double_t AliTPCFastMath::fgFastAsin[20000];
 158 AliTPCFastMath gAliTPCFastMath; // needed to fill the LUT
 159
 160 AliTPCFastMath::AliTPCFastMath(){
 161   //
 162   // initialized lookup table;
 163   for (Int_t i=0;i<10000;i++){
 164     fgFastAsin[2*i] = TMath::ASin(i/10000.);
 165     fgFastAsin[2*i+1] = (TMath::ASin((i+1)/10000.)-fgFastAsin[2*i]);
 166   }
 167 }
 168
 169 Double_t AliTPCFastMath::FastAsin(Double_t x){
 170   //
 171   // return asin using lookup table
 172   if (x>0){
 173     Int_t index = int(x*10000);
 174     return fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1];
 175   }
 176   x*=-1;
 177   Int_t index = int(x*10000);
 178   return -(fgFastAsin[2*index]+(x*10000.-index)*fgFastAsin[2*index+1]);
 179 }
 180 //__________________________________________________________________
 181 AliTPCtrackerMI::AliTPCtrackerMI()
 182                 :AliTracker(),
 183                  fkNIS(0),
 184                  fInnerSec(0),
 185                  fkNOS(0),
 186                  fOuterSec(0),
 187                  fN(0),
 188                  fSectors(0),
 189                  fInput(0),
 190                  fOutput(0),
 191                  fSeedTree(0),
 192                  fTreeDebug(0),
 193                  fEvent(0),
 194                  fDebug(0),
 195                  fNewIO(kFALSE),
 196                  fNtracks(0),
 197                  fSeeds(0),
 198                  fIteration(0),
 199                  fParam(0),
 200                  fDebugStreamer(0)
 201 {
 202   //
 203   // default constructor
 204   //
 205 }
 206 //_____________________________________________________________________
 207
 208
 209
 210 Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, Int_t accept){
 211   //
 212   //update track information using current cluster - track->fCurrentCluster
 213
 214
 215   AliTPCclusterMI* c =track->GetCurrentCluster();
 216   if (accept>0) track->SetCurrentClusterIndex1(track->GetCurrentClusterIndex1() | 0x8000);  //sign not accepted clusters
 217
 218   UInt_t i = track->GetCurrentClusterIndex1();
 219
 220   Int_t sec=(i&0xff000000)>>24;
 221   //Int_t row = (i&0x00ff0000)>>16;
 222   track->SetRow((i&0x00ff0000)>>16);
 223   track->SetSector(sec);
 224   //  Int_t index = i&0xFFFF;
 225   if (sec>=fParam->GetNInnerSector()) track->SetRow(track->GetRow()+fParam->GetNRowLow());
 226   track->SetClusterIndex2(track->GetRow(), i);
 227   //track->fFirstPoint = row;
 228   //if ( track->fLastPoint<row) track->fLastPoint =row;
 229   //  if (track->fRow<0 || track->fRow>160) {
 230   //  printf("problem\n");
 231   //}
 232   if (track->GetFirstPoint()>track->GetRow())
 233     track->SetFirstPoint(track->GetRow());
 234   if (track->GetLastPoint()<track->GetRow())
 235     track->SetLastPoint(track->GetRow());
 236
 237
 238   track->SetClusterPointer(track->GetRow(),c);
 239   //
 240
 241   Double_t angle2 = track->GetSnp()*track->GetSnp();
 242   //
 243   //SET NEW Track Point
 244   //
 245   if (angle2<1) //PH sometimes angle2 is very big. To be investigated...
 246   {
 247     angle2 = TMath::Sqrt(angle2/(1-angle2));
 248     AliTPCTrackerPoint   &point =*(track->GetTrackPoint(track->GetRow()));
 249     //
 250     point.SetSigmaY(c->GetSigmaY2()/track->GetCurrentSigmaY2());
 251     point.SetSigmaZ(c->GetSigmaZ2()/track->GetCurrentSigmaZ2());
 252     point.SetErrY(sqrt(track->GetErrorY2()));
 253     point.SetErrZ(sqrt(track->GetErrorZ2()));
 254     //
 255     point.SetX(track->GetX());
 256     point.SetY(track->GetY());
 257     point.SetZ(track->GetZ());
 258     point.SetAngleY(angle2);
 259     point.SetAngleZ(track->GetTgl());
 260     if (point.IsShared()){
 261       track->SetErrorY2(track->GetErrorY2()*4);
 262       track->SetErrorZ2(track->GetErrorZ2()*4);
 263     }
 264   }
 265
 266   Double_t chi2 = track->GetPredictedChi2(track->GetCurrentCluster());
 267   //
 268 //   track->SetErrorY2(track->GetErrorY2()*1.3);
 269 //   track->SetErrorY2(track->GetErrorY2()+0.01);
 270 //   track->SetErrorZ2(track->GetErrorZ2()*1.3);
 271 //   track->SetErrorZ2(track->GetErrorZ2()+0.005);
 272     //}
 273   if (accept>0) return 0;
 274   if (track->GetNumberOfClusters()%20==0){
 275     //    if (track->fHelixIn){
 276     //  TClonesArray & larr = *(track->fHelixIn);
 277     //  Int_t ihelix = larr.GetEntriesFast();
 278     //  new(larr[ihelix]) AliHelix(*track) ;
 279     //}
 280   }
 281   track->SetNoCluster(0);
 282   return track->Update(c,chi2,i);
 283 }
 284
 285
 286
 287 Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluster)
 288 {
 289   //
 290   // decide according desired precision to accept given
 291   // cluster for tracking
 292   Double_t sy2=ErrY2(seed,cluster);
 293   Double_t sz2=ErrZ2(seed,cluster);
 294
 295   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
 296   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
 297
 298   Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-seed->GetY())*
 299     (seed->GetCurrentCluster()->GetY()-seed->GetY())/sdistancey2;
 300   Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-seed->GetZ())*
 301     (seed->GetCurrentCluster()->GetZ()-seed->GetZ())/sdistancez2;
 302
 303   Double_t rdistance2  = rdistancey2+rdistancez2;
 304   //Int_t  accept =0;
 305
 306   if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
 307     Float_t rmsy2 = seed->GetCurrentSigmaY2();
 308     Float_t rmsz2 = seed->GetCurrentSigmaZ2();
 309     Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
 310     Float_t rmsz2p30 = seed->GetCMeanSigmaZ2p30();
 311     Float_t rmsy2p30R  = seed->GetCMeanSigmaY2p30R();
 312     Float_t rmsz2p30R  = seed->GetCMeanSigmaZ2p30R();
 313
 314     AliExternalTrackParam param(*seed);
 315     static TVectorD gcl(3),gtr(3);
 316     Float_t gclf[3];
 317     param.GetXYZ(gcl.GetMatrixArray());
 318     cluster->GetGlobalXYZ(gclf);
 319     gcl[0]=gclf[0];    gcl[1]=gclf[1];    gcl[2]=gclf[2];
 320
 321     if (AliTPCReconstructor::StreamLevel()>0) {
 322     (*fDebugStreamer)<<"ErrParam"<<
 323       "Cl.="<<cluster<<
 324       "T.="<<&param<<
 325       "gcl.="<<&gcl<<
 326       "gtr.="<<&gtr<<
 327       "erry2="<<sy2<<
 328       "errz2="<<sz2<<
 329       "rmsy2="<<rmsy2<<
 330       "rmsz2="<<rmsz2<<
 331       "rmsy2p30="<<rmsy2p30<<
 332       "rmsz2p30="<<rmsz2p30<<
 333       "rmsy2p30R="<<rmsy2p30R<<
 334       "rmsz2p30R="<<rmsz2p30R<<
 335       "\n";
 336     }
 337   }
 338
 339   if (rdistance2>16) return 3;
 340
 341
 342   if ((rdistancey2>9. || rdistancez2>9.) && cluster->GetType()==0)
 343     return 2;  //suspisiouce - will be changed
 344
 345   if ((rdistancey2>6.25 || rdistancez2>6.25) && cluster->GetType()>0)
 346     // strict cut on overlaped cluster
 347     return  2;  //suspisiouce - will be changed
 348
 349   if ( (rdistancey2>1. || rdistancez2>6.25 )
 350        && cluster->GetType()<0){
 351     seed->SetNFoundable(seed->GetNFoundable()-1);
 352     return 2;
 353   }
 354   return 0;
 355 }
 356
 357
 358
 359
 360 //_____________________________________________________________________________
 361 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par):
 362 AliTracker(),
 363                  fkNIS(par->GetNInnerSector()/2),
 364                  fInnerSec(0),
 365                  fkNOS(par->GetNOuterSector()/2),
 366                  fOuterSec(0),
 367                  fN(0),
 368                  fSectors(0),
 369                  fInput(0),
 370                  fOutput(0),
 371                  fSeedTree(0),
 372                  fTreeDebug(0),
 373                  fEvent(0),
 374                  fDebug(0),
 375                  fNewIO(0),
 376                  fNtracks(0),
 377                  fSeeds(0),
 378                  fIteration(0),
 379                  fParam(0),
 380                  fDebugStreamer(0)
 381 {
 382   //---------------------------------------------------------------------
 383   // The main TPC tracker constructor
 384   //---------------------------------------------------------------------
 385   fInnerSec=new AliTPCtrackerSector[fkNIS];
 386   fOuterSec=new AliTPCtrackerSector[fkNOS];
 387
 388   Int_t i;
 389   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
 390   for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
 391
 392   fParam = par;
 393   Int_t nrowlow = par->GetNRowLow();
 394   Int_t nrowup = par->GetNRowUp();
 395
 396
 397   for (Int_t i=0;i<nrowlow;i++){
 398     fXRow[i]     = par->GetPadRowRadiiLow(i);
 399     fPadLength[i]= par->GetPadPitchLength(0,i);
 400     fYMax[i]     = fXRow[i]*TMath::Tan(0.5*par->GetInnerAngle());
 401   }
 402
 403
 404   for (Int_t i=0;i<nrowup;i++){
 405     fXRow[i+nrowlow]      = par->GetPadRowRadiiUp(i);
 406     fPadLength[i+nrowlow] = par->GetPadPitchLength(60,i);
 407     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
 408   }
 409
 410   if (AliTPCReconstructor::StreamLevel()>0) {
 411     fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
 412   }
 413 }
 414 //________________________________________________________________________
 415 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
 416   AliTracker(t),
 417                  fkNIS(t.fkNIS),
 418                  fInnerSec(0),
 419                  fkNOS(t.fkNOS),
 420                  fOuterSec(0),
 421                  fN(0),
 422                  fSectors(0),
 423                  fInput(0),
 424                  fOutput(0),
 425                  fSeedTree(0),
 426                  fTreeDebug(0),
 427                  fEvent(0),
 428                  fDebug(0),
 429                  fNewIO(kFALSE),
 430                  fNtracks(0),
 431                  fSeeds(0),
 432                  fIteration(0),
 433                  fParam(0),
 434                  fDebugStreamer(0)
 435 {
 436   //------------------------------------
 437   // dummy copy constructor
 438   //------------------------------------------------------------------
 439   fOutput=t.fOutput;
 440 }
 441 AliTPCtrackerMI & AliTPCtrackerMI::operator=(const AliTPCtrackerMI& /*r*/){
 442   //------------------------------
 443   // dummy
 444   //--------------------------------------------------------------
 445   return *this;
 446 }
 447 //_____________________________________________________________________________
 448 AliTPCtrackerMI::~AliTPCtrackerMI() {
 449   //------------------------------------------------------------------
 450   // TPC tracker destructor
 451   //------------------------------------------------------------------
 452   delete[] fInnerSec;
 453   delete[] fOuterSec;
 454   if (fSeeds) {
 455     fSeeds->Delete();
 456     delete fSeeds;
 457   }
 458   if (fDebugStreamer) delete fDebugStreamer;
 459 }
 460
 461
 462 void AliTPCtrackerMI::FillESD(TObjArray* arr)
 463 {
 464   //
 465   //
 466   //fill esds using updated tracks
 467   if (fEvent){
 468     // write tracks to the event
 469     // store index of the track
 470     Int_t nseed=arr->GetEntriesFast();
 471     //FindKinks(arr,fEvent);
 472     for (Int_t i=0; i<nseed; i++) {
 473       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
 474       if (!pt) continue;
 475       pt->UpdatePoints();
 476       AddCovariance(pt);
 477       if (AliTPCReconstructor::StreamLevel()>1) {
 478         (*fDebugStreamer)<<"Track0"<<
 479           "Tr0.="<<pt<<
 480           "\n";
 481       }
 482       //      pt->PropagateTo(fParam->GetInnerRadiusLow());
 483       if (pt->GetKinkIndex(0)<=0){  //don't propagate daughter tracks
 484         pt->PropagateTo(fParam->GetInnerRadiusLow());
 485       }
 486
 487       if (( pt->GetPoints()[2]- pt->GetPoints()[0])>5 && pt->GetPoints()[3]>0.8){
 488         AliESDtrack iotrack;
 489         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 490         iotrack.SetTPCPoints(pt->GetPoints());
 491         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 492         iotrack.SetV0Indexes(pt->GetV0Indexes());
 493         //      iotrack.SetTPCpid(pt->fTPCr);
 494         //iotrack.SetTPCindex(i);
 495         fEvent->AddTrack(&iotrack);
 496         continue;
 497       }
 498
 499       if ( (pt->GetNumberOfClusters()>70)&& (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.55) {
 500         AliESDtrack iotrack;
 501         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 502         iotrack.SetTPCPoints(pt->GetPoints());
 503         //iotrack.SetTPCindex(i);
 504         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 505         iotrack.SetV0Indexes(pt->GetV0Indexes());
 506         //      iotrack.SetTPCpid(pt->fTPCr);
 507         fEvent->AddTrack(&iotrack);
 508         continue;
 509       }
 510       //
 511       // short tracks  - maybe decays
 512
 513       if ( (pt->GetNumberOfClusters()>30) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.70) {
 514         Int_t found,foundable,shared;
 515         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
 516         if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2)){
 517           AliESDtrack iotrack;
 518           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 519           //iotrack.SetTPCindex(i);
 520           iotrack.SetTPCPoints(pt->GetPoints());
 521           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 522           iotrack.SetV0Indexes(pt->GetV0Indexes());
 523           //iotrack.SetTPCpid(pt->fTPCr);
 524           fEvent->AddTrack(&iotrack);
 525           continue;
 526         }
 527       }
 528
 529       if ( (pt->GetNumberOfClusters()>20) && (Float_t(pt->GetNumberOfClusters())/Float_t(pt->GetNFoundable()))>0.8) {
 530         Int_t found,foundable,shared;
 531         pt->GetClusterStatistic(0,60,found, foundable,shared,kFALSE);
 532         if (found<20) continue;
 533         if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
 534         //
 535         AliESDtrack iotrack;
 536         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 537         iotrack.SetTPCPoints(pt->GetPoints());
 538         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 539         iotrack.SetV0Indexes(pt->GetV0Indexes());
 540         //iotrack.SetTPCpid(pt->fTPCr);
 541         //iotrack.SetTPCindex(i);
 542         fEvent->AddTrack(&iotrack);
 543         continue;
 544       }
 545       // short tracks  - secondaties
 546       //
 547       if ( (pt->GetNumberOfClusters()>30) ) {
 548         Int_t found,foundable,shared;
 549         pt->GetClusterStatistic(128,158,found, foundable,shared,kFALSE);
 550         if ( (found>20) && (pt->GetNShared()/float(pt->GetNumberOfClusters())<0.2) &&float(found)/float(foundable)>0.8){
 551           AliESDtrack iotrack;
 552           iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 553           iotrack.SetTPCPoints(pt->GetPoints());
 554           iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 555           iotrack.SetV0Indexes(pt->GetV0Indexes());
 556           //iotrack.SetTPCpid(pt->fTPCr);
 557           //iotrack.SetTPCindex(i);
 558           fEvent->AddTrack(&iotrack);
 559           continue;
 560         }
 561       }
 562
 563       if ( (pt->GetNumberOfClusters()>15)) {
 564         Int_t found,foundable,shared;
 565         pt->GetClusterStatistic(138,158,found, foundable,shared,kFALSE);
 566         if (found<15) continue;
 567         if (foundable<=0) continue;
 568         if (pt->GetNShared()/float(pt->GetNumberOfClusters())>0.2) continue;
 569         if (float(found)/float(foundable)<0.8) continue;
 570         //
 571         AliESDtrack iotrack;
 572         iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
 573         iotrack.SetTPCPoints(pt->GetPoints());
 574         iotrack.SetKinkIndexes(pt->GetKinkIndexes());
 575         iotrack.SetV0Indexes(pt->GetV0Indexes());
 576         //      iotrack.SetTPCpid(pt->fTPCr);
 577         //iotrack.SetTPCindex(i);
 578         fEvent->AddTrack(&iotrack);
 579         continue;
 580       }
 581     }
 582   }
 583   printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
 584 }
 585
 586
 587
 588
 589
 590 Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
 591   //
 592   //
 593   // Use calibrated cluster error from OCDB
 594   //
 595   AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
 596   //
 597   Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
 598   Int_t ctype = cl->GetType();
 599   Int_t    type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
 600   Double_t angle = seed->GetSnp()*seed->GetSnp();
 601   angle = TMath::Sqrt(TMath::Abs(angle/(1.-angle)));
 602   Double_t erry2 = clparam->GetError0Par(0,type, z,angle);
 603   if (ctype<0) {
 604     erry2+=0.5;  // edge cluster
 605   }
 606   erry2*=erry2;
 607   seed->SetErrorY2(erry2);
 608   //
 609   return erry2;
 610
 611 //calculate look-up table at the beginning
 612 //   static Bool_t  ginit = kFALSE;
 613 //   static Float_t gnoise1,gnoise2,gnoise3;
 614 //   static Float_t ggg1[10000];
 615 //   static Float_t ggg2[10000];
 616 //   static Float_t ggg3[10000];
 617 //   static Float_t glandau1[10000];
 618 //   static Float_t glandau2[10000];
 619 //   static Float_t glandau3[10000];
 620 //   //
 621 //   static Float_t gcor01[500];
 622 //   static Float_t gcor02[500];
 623 //   static Float_t gcorp[500];
 624 //   //
 625
 626 //   //
 627 //   if (ginit==kFALSE){
 628 //     for (Int_t i=1;i<500;i++){
 629 //       Float_t rsigma = float(i)/100.;
 630 //       gcor02[i] = TMath::Max(0.78 +TMath::Exp(7.4*(rsigma-1.2)),0.6);
 631 //       gcor01[i] = TMath::Max(0.72 +TMath::Exp(3.36*(rsigma-1.2)),0.6);
 632 //       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 633 //     }
 634
 635 //     //
 636 //     for (Int_t i=3;i<10000;i++){
 637 //       //
 638 //       //
 639 //       // inner sector
 640 //       Float_t amp = float(i);
 641 //       Float_t padlength =0.75;
 642 //       gnoise1 = 0.0004/padlength;
 643 //       Float_t nel     = 0.268*amp;
 644 //       Float_t nprim   = 0.155*amp;
 645 //       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 646 //       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 647 //       if (glandau1[i]>1) glandau1[i]=1;
 648 //       glandau1[i]*=padlength*padlength/12.;
 649 //       //
 650 //       // outer short
 651 //       padlength =1.;
 652 //       gnoise2   = 0.0004/padlength;
 653 //       nel       = 0.3*amp;
 654 //       nprim     = 0.133*amp;
 655 //       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 656 //       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 657 //       if (glandau2[i]>1) glandau2[i]=1;
 658 //       glandau2[i]*=padlength*padlength/12.;
 659 //       //
 660 //       //
 661 //       // outer long
 662 //       padlength =1.5;
 663 //       gnoise3   = 0.0004/padlength;
 664 //       nel       = 0.3*amp;
 665 //       nprim     = 0.133*amp;
 666 //       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 667 //       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 668 //       if (glandau3[i]>1) glandau3[i]=1;
 669 //       glandau3[i]*=padlength*padlength/12.;
 670 //       //
 671 //     }
 672 //     ginit = kTRUE;
 673 //   }
 674 //   //
 675 //   //
 676 //   //
 677 //   Int_t amp = int(TMath::Abs(cl->GetQ()));
 678 //   if (amp>9999) {
 679 //     seed->SetErrorY2(1.);
 680 //     return 1.;
 681 //   }
 682 //   Float_t snoise2;
 683 //   Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
 684 //   Int_t ctype = cl->GetType();
 685 //   Float_t padlength= GetPadPitchLength(seed->GetRow());
 686 //   Double_t angle2 = seed->GetSnp()*seed->GetSnp();
 687 //   angle2 = angle2/(1-angle2);
 688 //   //
 689 //   //cluster "quality"
 690 //   Int_t rsigmay = int(100.*cl->GetSigmaY2()/(seed->GetCurrentSigmaY2()));
 691 //   Float_t res;
 692 //   //
 693 //   if (fSectors==fInnerSec){
 694 //     snoise2 = gnoise1;
 695 //     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 696 //     if (ctype==0) res *= gcor01[rsigmay];
 697 //     if ((ctype>0)){
 698 //       res+=0.002;
 699 //       res*= gcorp[rsigmay];
 700 //     }
 701 //   }
 702 //   else {
 703 //     if (padlength<1.1){
 704 //       snoise2 = gnoise2;
 705 //       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 706 //       if (ctype==0) res *= gcor02[rsigmay];
 707 //       if ((ctype>0)){
 708 //      res+=0.002;
 709 //      res*= gcorp[rsigmay];
 710 //       }
 711 //     }
 712 //     else{
 713 //       snoise2 = gnoise3;
 714 //       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 715 //       if (ctype==0) res *= gcor02[rsigmay];
 716 //       if ((ctype>0)){
 717 //      res+=0.002;
 718 //      res*= gcorp[rsigmay];
 719 //       }
 720 //     }
 721 //   }
 722
 723 //   if (ctype<0){
 724 //     res+=0.005;
 725 //     res*=2.4;  // overestimate error 2 times
 726 //   }
 727 //   res+= snoise2;
 728
 729 //   if (res<2*snoise2)
 730 //     res = 2*snoise2;
 731
 732 //   seed->SetErrorY2(res);
 733 //   return res;
 734
 735
 736 }
 737
 738
 739
 740 Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 741   //
 742   //
 743   // Use calibrated cluster error from OCDB
 744   //
 745   AliTPCClusterParam * clparam = AliTPCcalibDB::Instance()->GetClusterParam();
 746   //
 747   Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
 748   Int_t ctype = cl->GetType();
 749   Int_t    type = (cl->GetRow()<63) ? 0: (cl->GetRow()>126) ? 1:2;
 750   //
 751   Double_t angle2 = seed->GetSnp()*seed->GetSnp();
 752   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
 753   Double_t angle = TMath::Sqrt(TMath::Abs(angle2));
 754   Double_t errz2 = clparam->GetError0Par(1,type, z,angle);
 755   if (ctype<0) {
 756     errz2+=0.5;  // edge cluster
 757   }
 758   errz2*=errz2;
 759   seed->SetErrorZ2(errz2);
 760   //
 761   return errz2;
 762
 763
 764
 765 //   //seed->SetErrorY2(0.1);
 766 //   //return 0.1;
 767 //   //calculate look-up table at the beginning
 768 //   static Bool_t  ginit = kFALSE;
 769 //   static Float_t gnoise1,gnoise2,gnoise3;
 770 //   static Float_t ggg1[10000];
 771 //   static Float_t ggg2[10000];
 772 //   static Float_t ggg3[10000];
 773 //   static Float_t glandau1[10000];
 774 //   static Float_t glandau2[10000];
 775 //   static Float_t glandau3[10000];
 776 //   //
 777 //   static Float_t gcor01[1000];
 778 //   static Float_t gcor02[1000];
 779 //   static Float_t gcorp[1000];
 780 //   //
 781
 782 //   //
 783 //   if (ginit==kFALSE){
 784 //     for (Int_t i=1;i<1000;i++){
 785 //       Float_t rsigma = float(i)/100.;
 786 //       gcor02[i] = TMath::Max(0.81 +TMath::Exp(6.8*(rsigma-1.2)),0.6);
 787 //       gcor01[i] = TMath::Max(0.72 +TMath::Exp(2.04*(rsigma-1.2)),0.6);
 788 //       gcorp[i]  = TMath::Max(TMath::Power((rsigma+0.5),1.5),1.2);
 789 //     }
 790
 791 //     //
 792 //     for (Int_t i=3;i<10000;i++){
 793 //       //
 794 //       //
 795 //       // inner sector
 796 //       Float_t amp = float(i);
 797 //       Float_t padlength =0.75;
 798 //       gnoise1 = 0.0004/padlength;
 799 //       Float_t nel     = 0.268*amp;
 800 //       Float_t nprim   = 0.155*amp;
 801 //       ggg1[i]          = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.001*nel/(padlength*padlength))/nel;
 802 //       glandau1[i]      = (2.+0.12*nprim)*0.5* (2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 803 //       if (glandau1[i]>1) glandau1[i]=1;
 804 //       glandau1[i]*=padlength*padlength/12.;
 805 //       //
 806 //       // outer short
 807 //       padlength =1.;
 808 //       gnoise2   = 0.0004/padlength;
 809 //       nel       = 0.3*amp;
 810 //       nprim     = 0.133*amp;
 811 //       ggg2[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 812 //       glandau2[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 813 //       if (glandau2[i]>1) glandau2[i]=1;
 814 //       glandau2[i]*=padlength*padlength/12.;
 815 //       //
 816 //       //
 817 //       // outer long
 818 //       padlength =1.5;
 819 //       gnoise3   = 0.0004/padlength;
 820 //       nel       = 0.3*amp;
 821 //       nprim     = 0.133*amp;
 822 //       ggg3[i]      = fParam->GetDiffT()*fParam->GetDiffT()*(2+0.0008*nel/(padlength*padlength))/nel;
 823 //       glandau3[i]  = (2.+0.12*nprim)*0.5*(2.+nprim*nprim*0.001/(padlength*padlength))/nprim;
 824 //       if (glandau3[i]>1) glandau3[i]=1;
 825 //       glandau3[i]*=padlength*padlength/12.;
 826 //       //
 827 //     }
 828 //     ginit = kTRUE;
 829 //   }
 830 //   //
 831 //   //
 832 //   //
 833 //   Int_t amp = int(TMath::Abs(cl->GetQ()));
 834 //   if (amp>9999) {
 835 //     seed->SetErrorY2(1.);
 836 //     return 1.;
 837 //   }
 838 //   Float_t snoise2;
 839 //   Float_t z = TMath::Abs(fParam->GetZLength(0)-TMath::Abs(seed->GetZ()));
 840 //   Int_t ctype = cl->GetType();
 841 //   Float_t padlength= GetPadPitchLength(seed->GetRow());
 842 //   //
 843 //   Double_t angle2 = seed->GetSnp()*seed->GetSnp();
 844 //   //  if (angle2<0.6) angle2 = 0.6;
 845 //   angle2 = seed->GetTgl()*seed->GetTgl()*(1+angle2/(1-angle2));
 846 //   //
 847 //   //cluster "quality"
 848 //   Int_t rsigmaz = int(100.*cl->GetSigmaZ2()/(seed->GetCurrentSigmaZ2()));
 849 //   Float_t res;
 850 //   //
 851 //   if (fSectors==fInnerSec){
 852 //     snoise2 = gnoise1;
 853 //     res     = ggg1[amp]*z+glandau1[amp]*angle2;
 854 //     if (ctype==0) res *= gcor01[rsigmaz];
 855 //     if ((ctype>0)){
 856 //       res+=0.002;
 857 //       res*= gcorp[rsigmaz];
 858 //     }
 859 //   }
 860 //   else {
 861 //     if (padlength<1.1){
 862 //       snoise2 = gnoise2;
 863 //       res     = ggg2[amp]*z+glandau2[amp]*angle2;
 864 //       if (ctype==0) res *= gcor02[rsigmaz];
 865 //       if ((ctype>0)){
 866 //      res+=0.002;
 867 //      res*= gcorp[rsigmaz];
 868 //       }
 869 //     }
 870 //     else{
 871 //       snoise2 = gnoise3;
 872 //       res     = ggg3[amp]*z+glandau3[amp]*angle2;
 873 //       if (ctype==0) res *= gcor02[rsigmaz];
 874 //       if ((ctype>0)){
 875 //      res+=0.002;
 876 //      res*= gcorp[rsigmaz];
 877 //       }
 878 //     }
 879 //   }
 880
 881 //   if (ctype<0){
 882 //     res+=0.002;
 883 //     res*=1.3;
 884 //   }
 885 //   if ((ctype<0) &&amp<70){
 886 //     res+=0.002;
 887 //     res*=1.3;
 888 //   }
 889 //   res += snoise2;
 890 //   if (res<2*snoise2)
 891 //      res = 2*snoise2;
 892 //   if (res>3) res =3;
 893 //   seed->SetErrorZ2(res);
 894 //   return res;
 895 }
 896
 897
 898
 899
 900
 901 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
 902 {
 903   //rotate to track "local coordinata
 904   Float_t x = seed->GetX();
 905   Float_t y = seed->GetY();
 906   Float_t ymax = x*TMath::Tan(0.5*fSectors->GetAlpha());
 907
 908   if (y > ymax) {
 909     seed->SetRelativeSector((seed->GetRelativeSector()+1) % fN);
 910     if (!seed->Rotate(fSectors->GetAlpha()))
 911       return;
 912   } else if (y <-ymax) {
 913     seed->SetRelativeSector((seed->GetRelativeSector()-1+fN) % fN);
 914     if (!seed->Rotate(-fSectors->GetAlpha()))
 915       return;
 916   }
 917
 918 }
 919
 920
 921
 922 //_____________________________________________________________________________
 923 Double_t AliTPCtrackerMI::F1old(Double_t x1,Double_t y1,
 924                    Double_t x2,Double_t y2,
 925                    Double_t x3,Double_t y3)
 926 {
 927   //-----------------------------------------------------------------
 928   // Initial approximation of the track curvature
 929   //-----------------------------------------------------------------
 930   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 931   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 932                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 933   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 934                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 935
 936   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 937   if ( xr*xr+yr*yr<=0.00000000000001) return 100;
 938   return -xr*yr/sqrt(xr*xr+yr*yr);
 939 }
 940
 941
 942
 943 //_____________________________________________________________________________
 944 Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
 945                    Double_t x2,Double_t y2,
 946                    Double_t x3,Double_t y3)
 947 {
 948   //-----------------------------------------------------------------
 949   // Initial approximation of the track curvature
 950   //-----------------------------------------------------------------
 951   x3 -=x1;
 952   x2 -=x1;
 953   y3 -=y1;
 954   y2 -=y1;
 955   //
 956   Double_t det = x3*y2-x2*y3;
 957   if (det==0) {
 958     return 100;
 959   }
 960   //
 961   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
 962   Double_t x0 = x3*0.5-y3*u;
 963   Double_t y0 = y3*0.5+x3*u;
 964   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
 965   if (det<0) c2*=-1;
 966   return c2;
 967 }
 968
 969
 970 Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
 971                    Double_t x2,Double_t y2,
 972                    Double_t x3,Double_t y3)
 973 {
 974   //-----------------------------------------------------------------
 975   // Initial approximation of the track curvature
 976   //-----------------------------------------------------------------
 977   x3 -=x1;
 978   x2 -=x1;
 979   y3 -=y1;
 980   y2 -=y1;
 981   //
 982   Double_t det = x3*y2-x2*y3;
 983   if (det==0) {
 984     return 100;
 985   }
 986   //
 987   Double_t u = 0.5* (x2*(x2-x3)+y2*(y2-y3))/det;
 988   Double_t x0 = x3*0.5-y3*u;
 989   Double_t y0 = y3*0.5+x3*u;
 990   Double_t c2 = 1/TMath::Sqrt(x0*x0+y0*y0);
 991   if (det<0) c2*=-1;
 992   x0+=x1;
 993   x0*=c2;
 994   return x0;
 995 }
 996
 997
 998
 999 //_____________________________________________________________________________
1000 Double_t AliTPCtrackerMI::F2old(Double_t x1,Double_t y1,
1001                    Double_t x2,Double_t y2,
1002                    Double_t x3,Double_t y3)
1003 {
1004   //-----------------------------------------------------------------
1005   // Initial approximation of the track curvature times center of curvature
1006   //-----------------------------------------------------------------
1007   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
1008   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
1009                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
1010   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
1011                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
1012
1013   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
1014
1015   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
1016 }
1017
1018 //_____________________________________________________________________________
1019 Double_t AliTPCtrackerMI::F3(Double_t x1,Double_t y1,
1020                    Double_t x2,Double_t y2,
1021                    Double_t z1,Double_t z2)
1022 {
1023   //-----------------------------------------------------------------
1024   // Initial approximation of the tangent of the track dip angle
1025   //-----------------------------------------------------------------
1026   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1027 }
1028
1029
1030 Double_t AliTPCtrackerMI::F3n(Double_t x1,Double_t y1,
1031                    Double_t x2,Double_t y2,
1032                    Double_t z1,Double_t z2, Double_t c)
1033 {
1034   //-----------------------------------------------------------------
1035   // Initial approximation of the tangent of the track dip angle
1036   //-----------------------------------------------------------------
1037
1038   //  Double_t angle1;
1039
1040   //angle1    =  (z1-z2)*c/(TMath::ASin(c*x1-ni)-TMath::ASin(c*x2-ni));
1041   //
1042   Double_t d  =  TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
1043   if (TMath::Abs(d*c*0.5)>1) return 0;
1044   //  Double_t   angle2    =  TMath::ASin(d*c*0.5);
1045   //  Double_t   angle2    =  AliTPCFastMath::FastAsin(d*c*0.5);
1046   Double_t   angle2    = (d*c*0.5>0.1)? TMath::ASin(d*c*0.5): AliTPCFastMath::FastAsin(d*c*0.5);
1047
1048   angle2  = (z1-z2)*c/(angle2*2.);
1049   return angle2;
1050 }
1051
1052 Bool_t   AliTPCtrackerMI::GetProlongation(Double_t x1, Double_t x2, Double_t x[5], Double_t &y, Double_t &z)
1053 {//-----------------------------------------------------------------
1054   // This function find proloncation of a track to a reference plane x=x2.
1055   //-----------------------------------------------------------------
1056
1057   Double_t dx=x2-x1;
1058
1059   if (TMath::Abs(x[4]*x1 - x[2]) >= 0.999) {
1060     return kFALSE;
1061   }
1062
1063   Double_t c1=x[4]*x1 - x[2], r1=sqrt(1.- c1*c1);
1064   Double_t c2=x[4]*x2 - x[2], r2=sqrt(1.- c2*c2);
1065   y = x[0];
1066   z = x[1];
1067
1068   Double_t dy = dx*(c1+c2)/(r1+r2);
1069   Double_t dz = 0;
1070   //
1071   Double_t delta = x[4]*dx*(c1+c2)/(c1*r2 + c2*r1);
1072
1073   if (TMath::Abs(delta)>0.01){
1074     dz = x[3]*TMath::ASin(delta)/x[4];
1075   }else{
1076     dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1077   }
1078
1079   //dz = x[3]*AliTPCFastMath::FastAsin(delta)/x[4];
1080
1081   y+=dy;
1082   z+=dz;
1083
1084   return kTRUE;
1085 }
1086
1087 Int_t  AliTPCtrackerMI::LoadClusters (TTree *tree)
1088 {
1089   //
1090   //
1091   fInput = tree;
1092   return LoadClusters();
1093 }
1094
1095
1096 Int_t  AliTPCtrackerMI::LoadClusters(TObjArray *arr)
1097 {
1098   //
1099   // load clusters to the memory
1100   AliTPCClustersRow *clrow = 0x0;
1101   Int_t lower   = arr->LowerBound();
1102   Int_t entries = arr->GetEntriesFast();
1103
1104   for (Int_t i=lower; i<entries; i++) {
1105     clrow = (AliTPCClustersRow*) arr->At(i);
1106         if(!clrow) continue;
1107         if(!clrow->GetArray()) continue;
1108
1109     //
1110     Int_t sec,row;
1111     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1112
1113     for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1114       Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1115     }
1116     //
1117     if (clrow->GetArray()->GetEntriesFast()<=0) continue;
1118     AliTPCtrackerRow * tpcrow=0;
1119     Int_t left=0;
1120     if (sec<fkNIS*2){
1121       tpcrow = &(fInnerSec[sec%fkNIS][row]);
1122       left = sec/fkNIS;
1123     }
1124     else{
1125       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1126       left = (sec-fkNIS*2)/fkNOS;
1127     }
1128     if (left ==0){
1129       tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1130       tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1131       for (Int_t i=0;i<tpcrow->GetN1();i++)
1132         tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1133     }
1134     if (left ==1){
1135       tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1136       tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1137       for (Int_t i=0;i<tpcrow->GetN2();i++)
1138         tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1139     }
1140   }
1141   //
1142   delete clrow;
1143   LoadOuterSectors();
1144   LoadInnerSectors();
1145   return 0;
1146 }
1147
1148 Int_t  AliTPCtrackerMI::LoadClusters(TClonesArray *arr)
1149 {
1150   //
1151   // load clusters to the memory from one
1152   // TClonesArray
1153   //
1154   AliTPCclusterMI *clust=0;
1155   Int_t count[72][96] = { {0} , {0} };
1156
1157   // loop over clusters
1158   for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1159     clust = (AliTPCclusterMI*)arr->At(icl);
1160     if(!clust) continue;
1161     //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
1162
1163     // transform clusters
1164     Transform(clust);
1165
1166     // count clusters per pad row
1167     count[clust->GetDetector()][clust->GetRow()]++;
1168   }
1169
1170   // insert clusters to sectors
1171   for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
1172     clust = (AliTPCclusterMI*)arr->At(icl);
1173     if(!clust) continue;
1174
1175     Int_t sec = clust->GetDetector();
1176     Int_t row = clust->GetRow();
1177
1178     // filter overlapping pad rows needed by HLT
1179     if(sec<fkNIS*2) { //IROCs
1180      if(row == 30) continue;
1181     }
1182     else { // OROCs
1183       if(row == 27 || row == 76) continue;
1184     }
1185
1186     Int_t left=0;
1187     if (sec<fkNIS*2){
1188       left = sec/fkNIS;
1189       fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fParam);
1190     }
1191     else{
1192       left = (sec-fkNIS*2)/fkNOS;
1193       fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fParam);
1194     }
1195   }
1196
1197   // Load functions must be called behind LoadCluster(TClonesArray*)
1198   // needed by HLT
1199   //LoadOuterSectors();
1200   //LoadInnerSectors();
1201
1202   return 0;
1203 }
1204
1205
1206 Int_t  AliTPCtrackerMI::LoadClusters()
1207 {
1208   //
1209   // load clusters to the memory
1210   AliTPCClustersRow *clrow= new AliTPCClustersRow;
1211   clrow->SetClass("AliTPCclusterMI");
1212   clrow->SetArray(0);
1213   clrow->GetArray()->ExpandCreateFast(10000);
1214   //
1215   //  TTree * tree = fClustersArray.GetTree();
1216
1217   TTree * tree = fInput;
1218   TBranch * br = tree->GetBranch("Segment");
1219   br->SetAddress(&clrow);
1220   //
1221   Int_t j=Int_t(tree->GetEntries());
1222   for (Int_t i=0; i<j; i++) {
1223     br->GetEntry(i);
1224     //
1225     Int_t sec,row;
1226     fParam->AdjustSectorRow(clrow->GetID(),sec,row);
1227     for (Int_t icl=0; icl<clrow->GetArray()->GetEntriesFast(); icl++){
1228       Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
1229     }
1230     //
1231     AliTPCtrackerRow * tpcrow=0;
1232     Int_t left=0;
1233     if (sec<fkNIS*2){
1234       tpcrow = &(fInnerSec[sec%fkNIS][row]);
1235       left = sec/fkNIS;
1236     }
1237     else{
1238       tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1239       left = (sec-fkNIS*2)/fkNOS;
1240     }
1241     if (left ==0){
1242       tpcrow->SetN1(clrow->GetArray()->GetEntriesFast());
1243       tpcrow->SetClusters1(new AliTPCclusterMI[tpcrow->GetN1()]);
1244       for (Int_t i=0;i<tpcrow->GetN1();i++)
1245         tpcrow->SetCluster1(i, *(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1246     }
1247     if (left ==1){
1248       tpcrow->SetN2(clrow->GetArray()->GetEntriesFast());
1249       tpcrow->SetClusters2(new AliTPCclusterMI[tpcrow->GetN2()]);
1250       for (Int_t i=0;i<tpcrow->GetN2();i++)
1251         tpcrow->SetCluster2(i,*(AliTPCclusterMI*)(clrow->GetArray()->At(i)));
1252     }
1253   }
1254   //
1255   delete clrow;
1256   LoadOuterSectors();
1257   LoadInnerSectors();
1258   return 0;
1259 }
1260
1261
1262 void AliTPCtrackerMI::UnloadClusters()
1263 {
1264   //
1265   // unload clusters from the memory
1266   //
1267   Int_t nrows = fOuterSec->GetNRows();
1268   for (Int_t sec = 0;sec<fkNOS;sec++)
1269     for (Int_t row = 0;row<nrows;row++){
1270       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1271       //      if (tpcrow){
1272       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1273       //        if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1274       //}
1275       tpcrow->ResetClusters();
1276     }
1277   //
1278   nrows = fInnerSec->GetNRows();
1279   for (Int_t sec = 0;sec<fkNIS;sec++)
1280     for (Int_t row = 0;row<nrows;row++){
1281       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1282       //if (tpcrow){
1283       //        if (tpcrow->fClusters1) delete []tpcrow->fClusters1;
1284       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2;
1285       //}
1286       tpcrow->ResetClusters();
1287     }
1288
1289   return ;
1290 }
1291
1292 void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
1293   //
1294   // Filling cluster to the array - For visualization purposes
1295   //
1296   Int_t nrows=0;
1297   nrows = fOuterSec->GetNRows();
1298   for (Int_t sec = 0;sec<fkNOS;sec++)
1299     for (Int_t row = 0;row<nrows;row++){
1300       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1301       if (!tpcrow) continue;
1302       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1303         array->AddLast((TObject*)((*tpcrow)[icl]));
1304       }
1305     }
1306   nrows = fInnerSec->GetNRows();
1307   for (Int_t sec = 0;sec<fkNIS;sec++)
1308     for (Int_t row = 0;row<nrows;row++){
1309       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1310       if (!tpcrow) continue;
1311       for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
1312         array->AddLast((TObject*)(*tpcrow)[icl]);
1313       }
1314     }
1315 }
1316
1317
1318 void   AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
1319   //
1320   //
1321   //
1322   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
1323   if (!transform) {
1324     AliFatal("Tranformations not in calibDB");
1325   }
1326   Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
1327   Int_t i[1]={cluster->GetDetector()};
1328   transform->Transform(x,i,0,1);
1329   if (cluster->GetDetector()%36>17){
1330     x[1]*=-1;
1331   }
1332
1333   //
1334   // in debug mode  check the transformation
1335   //
1336   if (AliTPCReconstructor::StreamLevel()>1) {
1337     Float_t gx[3];
1338     cluster->GetGlobalXYZ(gx);
1339     Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
1340     TTreeSRedirector &cstream = *fDebugStreamer;
1341     cstream<<"Transform"<<
1342       "event="<<event<<
1343       "x0="<<x[0]<<
1344       "x1="<<x[1]<<
1345       "x2="<<x[2]<<
1346       "gx0="<<gx[0]<<
1347       "gx1="<<gx[1]<<
1348       "gx2="<<gx[2]<<
1349       "Cl.="<<cluster<<
1350       "\n";
1351   }
1352   cluster->SetX(x[0]);
1353   cluster->SetY(x[1]);
1354   cluster->SetZ(x[2]);
1355   // The old stuff:
1356
1357   //
1358   //
1359   //
1360   //if (!fParam->IsGeoRead()) fParam->ReadGeoMatrices();
1361   TGeoHMatrix  *mat = fParam->GetClusterMatrix(cluster->GetDetector());
1362   //TGeoHMatrix  mat;
1363   Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1364   Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
1365   if (mat) mat->LocalToMaster(pos,posC);
1366   else{
1367     // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
1368   }
1369   cluster->SetX(posC[0]);
1370   cluster->SetY(posC[1]);
1371   cluster->SetZ(posC[2]);
1372 }
1373
1374 //_____________________________________________________________________________
1375 Int_t AliTPCtrackerMI::LoadOuterSectors() {
1376   //-----------------------------------------------------------------
1377   // This function fills outer TPC sectors with clusters.
1378   //-----------------------------------------------------------------
1379   Int_t nrows = fOuterSec->GetNRows();
1380   UInt_t index=0;
1381   for (Int_t sec = 0;sec<fkNOS;sec++)
1382     for (Int_t row = 0;row<nrows;row++){
1383       AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
1384       Int_t sec2 = sec+2*fkNIS;
1385       //left
1386       Int_t ncl = tpcrow->GetN1();
1387       while (ncl--) {
1388         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1389         index=(((sec2<<8)+row)<<16)+ncl;
1390         tpcrow->InsertCluster(c,index);
1391       }
1392       //right
1393       ncl = tpcrow->GetN2();
1394       while (ncl--) {
1395         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1396         index=((((sec2+fkNOS)<<8)+row)<<16)+ncl;
1397         tpcrow->InsertCluster(c,index);
1398       }
1399       //
1400       // write indexes for fast acces
1401       //
1402       for (Int_t i=0;i<510;i++)
1403         tpcrow->SetFastCluster(i,-1);
1404       for (Int_t i=0;i<tpcrow->GetN();i++){
1405         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1406         tpcrow->SetFastCluster(zi,i);  // write index
1407       }
1408       Int_t last = 0;
1409       for (Int_t i=0;i<510;i++){
1410         if (tpcrow->GetFastCluster(i)<0)
1411           tpcrow->SetFastCluster(i,last);
1412         else
1413           last = tpcrow->GetFastCluster(i);
1414       }
1415     }
1416   fN=fkNOS;
1417   fSectors=fOuterSec;
1418   return 0;
1419 }
1420
1421
1422 //_____________________________________________________________________________
1423 Int_t  AliTPCtrackerMI::LoadInnerSectors() {
1424   //-----------------------------------------------------------------
1425   // This function fills inner TPC sectors with clusters.
1426   //-----------------------------------------------------------------
1427   Int_t nrows = fInnerSec->GetNRows();
1428   UInt_t index=0;
1429   for (Int_t sec = 0;sec<fkNIS;sec++)
1430     for (Int_t row = 0;row<nrows;row++){
1431       AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
1432       //
1433       //left
1434       Int_t ncl = tpcrow->GetN1();
1435       while (ncl--) {
1436         AliTPCclusterMI *c= (tpcrow->GetCluster1(ncl));
1437         index=(((sec<<8)+row)<<16)+ncl;
1438         tpcrow->InsertCluster(c,index);
1439       }
1440       //right
1441       ncl = tpcrow->GetN2();
1442       while (ncl--) {
1443         AliTPCclusterMI *c= (tpcrow->GetCluster2(ncl));
1444         index=((((sec+fkNIS)<<8)+row)<<16)+ncl;
1445         tpcrow->InsertCluster(c,index);
1446       }
1447       //
1448       // write indexes for fast acces
1449       //
1450       for (Int_t i=0;i<510;i++)
1451         tpcrow->SetFastCluster(i,-1);
1452       for (Int_t i=0;i<tpcrow->GetN();i++){
1453         Int_t zi = Int_t((*tpcrow)[i]->GetZ()+255.);
1454         tpcrow->SetFastCluster(zi,i);  // write index
1455       }
1456       Int_t last = 0;
1457       for (Int_t i=0;i<510;i++){
1458         if (tpcrow->GetFastCluster(i)<0)
1459           tpcrow->SetFastCluster(i,last);
1460         else
1461           last = tpcrow->GetFastCluster(i);
1462       }
1463
1464     }
1465
1466   fN=fkNIS;
1467   fSectors=fInnerSec;
1468   return 0;
1469 }
1470
1471
1472
1473 //_________________________________________________________________________
1474 AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
1475   //--------------------------------------------------------------------
1476   //       Return pointer to a given cluster
1477   //--------------------------------------------------------------------
1478   if (index<0) return 0; // no cluster
1479   Int_t sec=(index&0xff000000)>>24;
1480   Int_t row=(index&0x00ff0000)>>16;
1481   Int_t ncl=(index&0x00007fff)>>00;
1482
1483   const AliTPCtrackerRow * tpcrow=0;
1484   AliTPCclusterMI * clrow =0;
1485
1486   if (sec<0 || sec>=fkNIS*4) {
1487     AliWarning(Form("Wrong sector %d",sec));
1488     return 0x0;
1489   }
1490
1491   if (sec<fkNIS*2){
1492     tpcrow = &(fInnerSec[sec%fkNIS][row]);
1493     if (tpcrow==0) return 0;
1494
1495     if (sec<fkNIS) {
1496       if (tpcrow->GetN1()<=ncl) return 0;
1497       clrow = tpcrow->GetClusters1();
1498     }
1499     else {
1500       if (tpcrow->GetN2()<=ncl) return 0;
1501       clrow = tpcrow->GetClusters2();
1502     }
1503   }
1504   else {
1505     tpcrow = &(fOuterSec[(sec-fkNIS*2)%fkNOS][row]);
1506     if (tpcrow==0) return 0;
1507
1508     if (sec-2*fkNIS<fkNOS) {
1509       if (tpcrow->GetN1()<=ncl) return 0;
1510       clrow = tpcrow->GetClusters1();
1511     }
1512     else {
1513       if (tpcrow->GetN2()<=ncl) return 0;
1514       clrow = tpcrow->GetClusters2();
1515     }
1516   }
1517
1518   return &(clrow[ncl]);
1519
1520 }
1521
1522
1523
1524 Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
1525   //-----------------------------------------------------------------
1526   // This function tries to find a track prolongation to next pad row
1527   //-----------------------------------------------------------------
1528   //
1529   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1530   AliTPCclusterMI *cl=0;
1531   Int_t tpcindex= t.GetClusterIndex2(nr);
1532   //
1533   // update current shape info every 5 pad-row
1534   //  if ( (nr%5==0) || t.GetNumberOfClusters()<2 || (t.fCurrentSigmaY2<0.0001) ){
1535     GetShape(&t,nr);
1536     //}
1537   //
1538   if (fIteration>0 && tpcindex>=-1){  //if we have already clusters
1539     //
1540     if (tpcindex==-1) return 0; //track in dead zone
1541     if (tpcindex>0){     //
1542       cl = t.GetClusterPointer(nr);
1543       if ( (cl==0) ) cl = GetClusterMI(tpcindex);
1544       t.SetCurrentClusterIndex1(tpcindex);
1545     }
1546     if (cl){
1547       Int_t relativesector = ((tpcindex&0xff000000)>>24)%18;  // if previously accepted cluster in different sector
1548       Float_t angle = relativesector*fSectors->GetAlpha()+fSectors->GetAlphaShift();
1549       //
1550       if (angle<-TMath::Pi()) angle += 2*TMath::Pi();
1551       if (angle>=TMath::Pi()) angle -= 2*TMath::Pi();
1552
1553       if (TMath::Abs(angle-t.GetAlpha())>0.001){
1554         Double_t rotation = angle-t.GetAlpha();
1555         t.SetRelativeSector(relativesector);
1556         if (!t.Rotate(rotation)) return 0;
1557       }
1558       if (!t.PropagateTo(x)) return 0;
1559       //
1560       t.SetCurrentCluster(cl);
1561       t.SetRow(nr);
1562       Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1563       if ((tpcindex&0x8000)==0) accept =0;
1564       if (accept<3) {
1565         //if founded cluster is acceptible
1566         if (cl->IsUsed(11)) {  // id cluster is shared inrease uncertainty
1567           t.SetErrorY2(t.GetErrorY2()+0.03);
1568           t.SetErrorZ2(t.GetErrorZ2()+0.03);
1569           t.SetErrorY2(t.GetErrorY2()*3);
1570           t.SetErrorZ2(t.GetErrorZ2()*3);
1571         }
1572         t.SetNFoundable(t.GetNFoundable()+1);
1573         UpdateTrack(&t,accept);
1574         return 1;
1575       }
1576     }
1577   }
1578   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;  // cut on angle
1579   if (fIteration>1){
1580     // not look for new cluster during refitting
1581     t.SetNFoundable(t.GetNFoundable()+1);
1582     return 0;
1583   }
1584   //
1585   UInt_t index=0;
1586   //  if (TMath::Abs(t.GetSnp())>0.95 || TMath::Abs(x*t.GetC()-t.GetEta())>0.95) return 0;// patch 28 fev 06
1587   Double_t  y=t.GetYat(x);
1588   if (TMath::Abs(y)>ymax){
1589     if (y > ymax) {
1590       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1591       if (!t.Rotate(fSectors->GetAlpha()))
1592         return 0;
1593     } else if (y <-ymax) {
1594       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1595       if (!t.Rotate(-fSectors->GetAlpha()))
1596         return 0;
1597     }
1598     //return 1;
1599   }
1600   //
1601   if (!t.PropagateTo(x)) {
1602     if (fIteration==0) t.SetRemoval(10);
1603     return 0;
1604   }
1605   y=t.GetY();
1606   Double_t z=t.GetZ();
1607   //
1608
1609   if (!IsActive(t.GetRelativeSector(),nr)) {
1610     t.SetInDead(kTRUE);
1611     t.SetClusterIndex2(nr,-1);
1612     return 0;
1613   }
1614   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1615   Bool_t isActive  = IsActive(t.GetRelativeSector(),nr);
1616   Bool_t isActive2 = (nr>=fInnerSec->GetNRows()) ? fOuterSec[t.GetRelativeSector()][nr-fInnerSec->GetNRows()].GetN()>0:fInnerSec[t.GetRelativeSector()][nr].GetN()>0;
1617
1618   if (!isActive || !isActive2) return 0;
1619
1620   const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1621   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1622   Double_t  roady  =1.;
1623   Double_t  roadz = 1.;
1624   //
1625   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1626     t.SetInDead(kTRUE);
1627     t.SetClusterIndex2(nr,-1);
1628     return 0;
1629   }
1630   else
1631     {
1632       if (TMath::Abs(z)<(AliTPCReconstructor::GetCtgRange()*x+10) && TMath::Abs(z)<fParam->GetZLength(0) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker()))
1633         t.SetNFoundable(t.GetNFoundable()+1);
1634       else
1635         return 0;
1636     }
1637   //calculate
1638   if (krow) {
1639     //    cl = krow.FindNearest2(y+10.,z,roady,roadz,index);
1640     cl = krow.FindNearest2(y,z,roady,roadz,index);
1641     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1642   }
1643   if (cl) {
1644     t.SetCurrentCluster(cl);
1645     t.SetRow(nr);
1646     if (fIteration==2&&cl->IsUsed(10)) return 0;
1647     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1648     if (fIteration==2&&cl->IsUsed(11)) {
1649       t.SetErrorY2(t.GetErrorY2()+0.03);
1650       t.SetErrorZ2(t.GetErrorZ2()+0.03);
1651       t.SetErrorY2(t.GetErrorY2()*3);
1652       t.SetErrorZ2(t.GetErrorZ2()*3);
1653     }
1654     /*
1655     if (t.fCurrentCluster->IsUsed(10)){
1656       //
1657       //
1658
1659       t.fNShared++;
1660       if (t.fNShared>0.7*t.GetNumberOfClusters()) {
1661         t.fRemoval =10;
1662         return 0;
1663       }
1664     }
1665     */
1666     if (accept<3) UpdateTrack(&t,accept);
1667
1668   } else {
1669     if ( fIteration==0 && t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) t.SetRemoval(10);
1670
1671   }
1672   return 1;
1673 }
1674
1675 Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
1676   //-----------------------------------------------------------------
1677   // This function tries to find a track prolongation to next pad row
1678   //-----------------------------------------------------------------
1679   //
1680   Double_t  x= GetXrow(nr), ymax=GetMaxY(nr);
1681   Double_t y,z;
1682   if (!t.GetProlongation(x,y,z)) {
1683     t.SetRemoval(10);
1684     return 0;
1685   }
1686   //
1687   //
1688   if (TMath::Abs(y)>ymax){
1689
1690     if (y > ymax) {
1691       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1692       if (!t.Rotate(fSectors->GetAlpha()))
1693         return 0;
1694     } else if (y <-ymax) {
1695       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1696       if (!t.Rotate(-fSectors->GetAlpha()))
1697         return 0;
1698     }
1699     if (!t.PropagateTo(x)) {
1700       return 0;
1701     }
1702     t.GetProlongation(x,y,z);
1703   }
1704   //
1705   // update current shape info every 2 pad-row
1706   if ( (nr%2==0) || t.GetNumberOfClusters()<2 || (t.GetCurrentSigmaY2()<0.0001) ){
1707     //    t.fCurrentSigmaY = GetSigmaY(&t);
1708     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1709     GetShape(&t,nr);
1710   }
1711   //
1712   AliTPCclusterMI *cl=0;
1713   UInt_t index=0;
1714
1715
1716   //Int_t nr2 = nr;
1717   const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1718   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
1719   Double_t  roady  =1.;
1720   Double_t  roadz = 1.;
1721   //
1722   Int_t row = nr;
1723   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1724     t.SetInDead(kTRUE);
1725     t.SetClusterIndex2(row,-1);
1726     return 0;
1727   }
1728   else
1729     {
1730       if (TMath::Abs(z)>(AliTPCReconstructor::GetCtgRange()*x+10)) t.SetClusterIndex2(row,-1);
1731     }
1732   //calculate
1733
1734   if ((cl==0)&&(krow)) {
1735     //    cl = krow.FindNearest2(y+10,z,roady,roadz,index);
1736     cl = krow.FindNearest2(y,z,roady,roadz,index);
1737
1738     if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(index));
1739   }
1740
1741   if (cl) {
1742     t.SetCurrentCluster(cl);
1743     //    Int_t accept = AcceptCluster(&t,t.fCurrentCluster);
1744     //if (accept<3){
1745       t.SetClusterIndex2(row,index);
1746       t.SetClusterPointer(row, cl);
1747       //}
1748   }
1749   return 1;
1750 }
1751
1752
1753
1754 //_________________________________________________________________________
1755 Bool_t AliTPCtrackerMI::GetTrackPoint(Int_t index, AliTrackPoint &p ) const
1756 {
1757   // Get track space point by index
1758   // return false in case the cluster doesn't exist
1759   AliTPCclusterMI *cl = GetClusterMI(index);
1760   if (!cl) return kFALSE;
1761   Int_t sector = (index&0xff000000)>>24;
1762   //  Int_t row = (index&0x00ff0000)>>16;
1763   Float_t xyz[3];
1764   //  xyz[0] = fParam->GetPadRowRadii(sector,row);
1765   xyz[0] = cl->GetX();
1766   xyz[1] = cl->GetY();
1767   xyz[2] = cl->GetZ();
1768   Float_t sin,cos;
1769   fParam->AdjustCosSin(sector,cos,sin);
1770   Float_t x = cos*xyz[0]-sin*xyz[1];
1771   Float_t y = cos*xyz[1]+sin*xyz[0];
1772   Float_t cov[6];
1773   Float_t sigmaY2 = 0.027*cl->GetSigmaY2();
1774   if (sector < fParam->GetNInnerSector()) sigmaY2 *= 2.07;
1775   Float_t sigmaZ2 = 0.066*cl->GetSigmaZ2();
1776   if (sector < fParam->GetNInnerSector()) sigmaZ2 *= 1.77;
1777   cov[0] = sin*sin*sigmaY2;
1778   cov[1] = -sin*cos*sigmaY2;
1779   cov[2] = 0.;
1780   cov[3] = cos*cos*sigmaY2;
1781   cov[4] = 0.;
1782   cov[5] = sigmaZ2;
1783   p.SetXYZ(x,y,xyz[2],cov);
1784   AliGeomManager::ELayerID iLayer;
1785   Int_t idet;
1786   if (sector < fParam->GetNInnerSector()) {
1787     iLayer = AliGeomManager::kTPC1;
1788     idet = sector;
1789   }
1790   else {
1791     iLayer = AliGeomManager::kTPC2;
1792     idet = sector - fParam->GetNInnerSector();
1793   }
1794   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
1795   p.SetVolumeID(volid);
1796   return kTRUE;
1797 }
1798
1799
1800
1801 Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
1802   //-----------------------------------------------------------------
1803   // This function tries to find a track prolongation to next pad row
1804   //-----------------------------------------------------------------
1805   t.SetCurrentCluster(0);
1806   t.SetCurrentClusterIndex1(0);
1807
1808   Double_t xt=t.GetX();
1809   Int_t     row = GetRowNumber(xt)-1;
1810   Double_t  ymax= GetMaxY(nr);
1811
1812   if (row < nr) return 1; // don't prolongate if not information until now -
1813 //   if (TMath::Abs(t.GetSnp())>0.9 && t.GetNumberOfClusters()>40. && fIteration!=2) {
1814 //     t.fRemoval =10;
1815 //     return 0;  // not prolongate strongly inclined tracks
1816 //   }
1817 //   if (TMath::Abs(t.GetSnp())>0.95) {
1818 //     t.fRemoval =10;
1819 //     return 0;  // not prolongate strongly inclined tracks
1820 //   }// patch 28 fev 06
1821
1822   Double_t x= GetXrow(nr);
1823   Double_t y,z;
1824   //t.PropagateTo(x+0.02);
1825   //t.PropagateTo(x+0.01);
1826   if (!t.PropagateTo(x)){
1827     return 0;
1828   }
1829   //
1830   y=t.GetY();
1831   z=t.GetZ();
1832
1833   if (TMath::Abs(y)>ymax){
1834     if (y > ymax) {
1835       t.SetRelativeSector((t.GetRelativeSector()+1) % fN);
1836       if (!t.Rotate(fSectors->GetAlpha()))
1837         return 0;
1838     } else if (y <-ymax) {
1839       t.SetRelativeSector((t.GetRelativeSector()-1+fN) % fN);
1840       if (!t.Rotate(-fSectors->GetAlpha()))
1841         return 0;
1842     }
1843     //    if (!t.PropagateTo(x)){
1844     //  return 0;
1845     //}
1846     return 1;
1847     //y = t.GetY();
1848   }
1849   //
1850   if (TMath::Abs(t.GetSnp())>AliTPCReconstructor::GetMaxSnpTracker()) return 0;
1851
1852   if (!IsActive(t.GetRelativeSector(),nr)) {
1853     t.SetInDead(kTRUE);
1854     t.SetClusterIndex2(nr,-1);
1855     return 0;
1856   }
1857   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
1858
1859   AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
1860
1861   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
1862     t.SetInDead(kTRUE);
1863     t.SetClusterIndex2(nr,-1);
1864     return 0;
1865   }
1866   else
1867     {
1868       if (TMath::Abs(t.GetZ())<(AliTPCReconstructor::GetCtgRange()*t.GetX()+10) && (TMath::Abs(t.GetSnp())<AliTPCReconstructor::GetMaxSnpTracker())) t.SetNFoundable(t.GetNFoundable()+1);
1869       else
1870         return 0;
1871     }
1872
1873   // update current
1874   if ( (nr%2==0) || t.GetNumberOfClusters()<2){
1875     //    t.fCurrentSigmaY = GetSigmaY(&t);
1876     //t.fCurrentSigmaZ = GetSigmaZ(&t);
1877     GetShape(&t,nr);
1878   }
1879
1880   AliTPCclusterMI *cl=0;
1881   Int_t index=0;
1882   //
1883   Double_t roady = 1.;
1884   Double_t roadz = 1.;
1885   //
1886
1887   if (!cl){
1888     index = t.GetClusterIndex2(nr);
1889     if ( (index>0) && (index&0x8000)==0){
1890       cl = t.GetClusterPointer(nr);
1891       if ( (cl==0) && (index>0)) cl = GetClusterMI(index);
1892       t.SetCurrentClusterIndex1(index);
1893       if (cl) {
1894         t.SetCurrentCluster(cl);
1895         return 1;
1896       }
1897     }
1898   }
1899
1900   //  if (index<0) return 0;
1901   UInt_t uindex = TMath::Abs(index);
1902
1903   if (krow) {
1904     //cl = krow.FindNearest2(y+10,z,roady,roadz,uindex);
1905     cl = krow.FindNearest2(y,z,roady,roadz,uindex);
1906   }
1907
1908   if (cl) t.SetCurrentClusterIndex1(krow.GetIndex(uindex));
1909   t.SetCurrentCluster(cl);
1910
1911   return 1;
1912 }
1913
1914
1915 Int_t AliTPCtrackerMI::FollowToNextCluster(AliTPCseed & t, Int_t nr) {
1916   //-----------------------------------------------------------------
1917   // This function tries to find a track prolongation to next pad row
1918   //-----------------------------------------------------------------
1919
1920   //update error according neighborhoud
1921
1922   if (t.GetCurrentCluster()) {
1923     t.SetRow(nr);
1924     Int_t accept = AcceptCluster(&t,t.GetCurrentCluster());
1925
1926     if (t.GetCurrentCluster()->IsUsed(10)){
1927       //
1928       //
1929       //  t.fErrorZ2*=2;
1930       //  t.fErrorY2*=2;
1931       t.SetNShared(t.GetNShared()+1);
1932       if (t.GetNShared()>0.7*t.GetNumberOfClusters()) {
1933         t.SetRemoval(10);
1934         return 0;
1935       }
1936     }
1937     if (fIteration>0) accept = 0;
1938     if (accept<3)  UpdateTrack(&t,accept);
1939
1940   } else {
1941     if (fIteration==0){
1942       if ( ( (t.GetSigmaY2()+t.GetSigmaZ2())>0.16)&& t.GetNumberOfClusters()>18) t.SetRemoval(10);
1943       if (  t.GetChi2()/t.GetNumberOfClusters()>6 &&t.GetNumberOfClusters()>18) t.SetRemoval(10);
1944
1945       if (( (t.GetNFoundable()*0.5 > t.GetNumberOfClusters()) || t.GetNoCluster()>15)) t.SetRemoval(10);
1946     }
1947   }
1948   return 1;
1949 }
1950
1951
1952
1953 //_____________________________________________________________________________
1954 Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf, Int_t step) {
1955   //-----------------------------------------------------------------
1956   // This function tries to find a track prolongation.
1957   //-----------------------------------------------------------------
1958   Double_t xt=t.GetX();
1959   //
1960   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
1961   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
1962   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
1963   //
1964   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
1965
1966   Int_t first = GetRowNumber(xt)-1;
1967   for (Int_t nr= first; nr>=rf; nr-=step) {
1968     // update kink info
1969     if (t.GetKinkIndexes()[0]>0){
1970       for (Int_t i=0;i<3;i++){
1971         Int_t index = t.GetKinkIndexes()[i];
1972         if (index==0) break;
1973         if (index<0) continue;
1974         //
1975         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
1976         if (!kink){
1977           printf("PROBLEM\n");
1978         }
1979         else{
1980           Int_t kinkrow = kink->GetTPCRow0()+2+Int_t(0.5/(0.05+kink->GetAngle(2)));
1981           if (kinkrow==nr){
1982             AliExternalTrackParam paramd(t);
1983             kink->SetDaughter(paramd);
1984             kink->SetStatus(2,5);
1985             kink->Update();
1986           }
1987         }
1988       }
1989     }
1990
1991     if (nr==80) t.UpdateReference();
1992     if (nr<fInnerSec->GetNRows())
1993       fSectors = fInnerSec;
1994     else
1995       fSectors = fOuterSec;
1996     if (FollowToNext(t,nr)==0)
1997       if (!t.IsActive())
1998         return 0;
1999
2000   }
2001   return 1;
2002 }
2003
2004
2005 //_____________________________________________________________________________
2006 Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t rf, Int_t step) {
2007   //-----------------------------------------------------------------
2008   // This function tries to find a track prolongation.
2009   //-----------------------------------------------------------------
2010   Double_t xt=t.GetX();
2011   //
2012   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2013   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2014   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2015   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2016
2017   for (Int_t nr=GetRowNumber(xt)-1; nr>=rf; nr-=step) {
2018
2019     if (FollowToNextFast(t,nr)==0)
2020       if (!t.IsActive()) return 0;
2021
2022   }
2023   return 1;
2024 }
2025
2026
2027
2028
2029
2030 Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
2031   //-----------------------------------------------------------------
2032   // This function tries to find a track prolongation.
2033   //-----------------------------------------------------------------
2034   //
2035   Double_t xt=t.GetX();
2036   Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
2037   if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2038   if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2039   t.SetRelativeSector(Int_t(alpha/fSectors->GetAlpha()+0.0001)%fN);
2040
2041   Int_t first = t.GetFirstPoint();
2042   if (first<GetRowNumber(xt)+1) first = GetRowNumber(xt)+1;
2043   //
2044   if (first<0) first=0;
2045   for (Int_t nr=first; nr<=rf; nr++) {
2046     //    if ( (TMath::Abs(t.GetSnp())>0.95)) break;//patch 28 fev 06
2047     if (t.GetKinkIndexes()[0]<0){
2048       for (Int_t i=0;i<3;i++){
2049         Int_t index = t.GetKinkIndexes()[i];
2050         if (index==0) break;
2051         if (index>0) continue;
2052         index = TMath::Abs(index);
2053         AliKink * kink = (AliKink*)fEvent->GetKink(index-1);
2054         if (!kink){
2055           printf("PROBLEM\n");
2056         }
2057         else{
2058           Int_t kinkrow = kink->GetTPCRow0()-2-Int_t(0.5/(0.05+kink->GetAngle(2)));
2059           if (kinkrow==nr){
2060             AliExternalTrackParam paramm(t);
2061             kink->SetMother(paramm);
2062             kink->SetStatus(2,1);
2063             kink->Update();
2064           }
2065         }
2066       }
2067     }
2068     //
2069     if (nr<fInnerSec->GetNRows())
2070       fSectors = fInnerSec;
2071     else
2072       fSectors = fOuterSec;
2073
2074     FollowToNext(t,nr);
2075   }
2076   return 1;
2077 }
2078
2079
2080
2081
2082
2083 Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
2084 {
2085   //
2086   //
2087   sum1=0;
2088   sum2=0;
2089   Int_t sum=0;
2090   //
2091   Float_t dz2 =(s1->GetZ() - s2->GetZ());
2092   dz2*=dz2;
2093
2094   Float_t dy2 =TMath::Abs((s1->GetY() - s2->GetY()));
2095   dy2*=dy2;
2096   Float_t distance = TMath::Sqrt(dz2+dy2);
2097   if (distance>4.) return 0; // if there are far away  - not overlap - to reduce combinatorics
2098
2099   //  Int_t offset =0;
2100   Int_t firstpoint = TMath::Min(s1->GetFirstPoint(),s2->GetFirstPoint());
2101   Int_t lastpoint = TMath::Max(s1->GetLastPoint(),s2->GetLastPoint());
2102   if (lastpoint>160)
2103     lastpoint =160;
2104   if (firstpoint<0)
2105     firstpoint = 0;
2106   if (firstpoint>lastpoint) {
2107     firstpoint =lastpoint;
2108     //    lastpoint  =160;
2109   }
2110
2111
2112   for (Int_t i=firstpoint-1;i<lastpoint+1;i++){
2113     if (s1->GetClusterIndex2(i)>0) sum1++;
2114     if (s2->GetClusterIndex2(i)>0) sum2++;
2115     if (s1->GetClusterIndex2(i)==s2->GetClusterIndex2(i) && s1->GetClusterIndex2(i)>0) {
2116       sum++;
2117     }
2118   }
2119   if (sum<5) return 0;
2120
2121   Float_t summin = TMath::Min(sum1+1,sum2+1);
2122   Float_t ratio = (sum+1)/Float_t(summin);
2123   return ratio;
2124 }
2125
2126 void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
2127 {
2128   //
2129   //
2130   Float_t thetaCut = 0.2;//+10.*TMath::Sqrt(s1->GetSigmaTglZ()+ s2->GetSigmaTglZ());
2131   if (TMath::Abs(s1->GetTgl()-s2->GetTgl())>thetaCut) return;
2132   Float_t minCl = TMath::Min(s1->GetNumberOfClusters(),s2->GetNumberOfClusters());
2133   Int_t cutN0 = TMath::Max(5,TMath::Nint(0.1*minCl));
2134
2135   //
2136   Int_t sumshared=0;
2137   //
2138   //Int_t firstpoint = TMath::Max(s1->GetFirstPoint(),s2->GetFirstPoint());
2139   //Int_t lastpoint = TMath::Min(s1->GetLastPoint(),s2->GetLastPoint());
2140   Int_t firstpoint = 0;
2141   Int_t lastpoint = 160;
2142   //
2143   //  if (firstpoint>=lastpoint-5) return;;
2144
2145   for (Int_t i=firstpoint;i<lastpoint;i++){
2146     //    if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2147     if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2148       sumshared++;
2149       s1->SetSharedMapBit(i, kTRUE);
2150       s2->SetSharedMapBit(i, kTRUE);
2151     }
2152     if (s1->GetClusterIndex2(i)>0)
2153       s1->SetClusterMapBit(i, kTRUE);
2154   }
2155   if (sumshared>cutN0){
2156     // sign clusters
2157     //
2158     for (Int_t i=firstpoint;i<lastpoint;i++){
2159       //      if ( (s1->GetClusterIndex2(i)&0xFFFF8FFF)==(s2->GetClusterIndex2(i)&0xFFFF8FFF) && s1->GetClusterIndex2(i)>0) {
2160       if ( (s1->GetClusterIndex2(i))==(s2->GetClusterIndex2(i)) && s1->GetClusterIndex2(i)>0) {
2161         AliTPCTrackerPoint *p1  = s1->GetTrackPoint(i);
2162         AliTPCTrackerPoint *p2  = s2->GetTrackPoint(i);;
2163         if (s1->IsActive()&&s2->IsActive()){
2164           p1->SetShared(kTRUE);
2165           p2->SetShared(kTRUE);
2166         }
2167       }
2168     }
2169   }
2170   //
2171   if (sumshared>cutN0){
2172     for (Int_t i=0;i<4;i++){
2173       if (s1->GetOverlapLabel(3*i)==0){
2174         s1->SetOverlapLabel(3*i,  s2->GetLabel());
2175         s1->SetOverlapLabel(3*i+1,sumshared);
2176         s1->SetOverlapLabel(3*i+2,s2->GetUniqueID());
2177         break;
2178       }
2179     }
2180     for (Int_t i=0;i<4;i++){
2181       if (s2->GetOverlapLabel(3*i)==0){
2182         s2->SetOverlapLabel(3*i,  s1->GetLabel());
2183         s2->SetOverlapLabel(3*i+1,sumshared);
2184         s2->SetOverlapLabel(3*i+2,s1->GetUniqueID());
2185         break;
2186       }
2187     }
2188   }
2189 }
2190
2191 void  AliTPCtrackerMI::SignShared(TObjArray * arr)
2192 {
2193   //
2194   //sort trackss according sectors
2195   //
2196   for (Int_t i=0; i<arr->GetEntriesFast(); i++) {
2197     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2198     if (!pt) continue;
2199     //if (pt) RotateToLocal(pt);
2200     pt->SetSort(0);
2201   }
2202   arr->UnSort();
2203   arr->Sort();  // sorting according relative sectors
2204   arr->Expand(arr->GetEntries());
2205   //
2206   //
2207   Int_t nseed=arr->GetEntriesFast();
2208   for (Int_t i=0; i<nseed; i++) {
2209     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2210     if (!pt) continue;
2211     for (Int_t j=0;j<=12;j++){
2212       pt->SetOverlapLabel(j,0);
2213     }
2214   }
2215   for (Int_t i=0; i<nseed; i++) {
2216     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2217     if (!pt) continue;
2218     if (pt->GetRemoval()>10) continue;
2219     for (Int_t j=i+1; j<nseed; j++){
2220       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
2221       if (TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>1) continue;
2222       //      if (pt2){
2223       if (pt2->GetRemoval()<=10) {
2224         //if ( TMath::Abs(pt->GetRelativeSector()-pt2->GetRelativeSector())>0) break;
2225         SignShared(pt,pt2);
2226       }
2227     }
2228   }
2229 }
2230
2231
2232 void AliTPCtrackerMI::SortTracks(TObjArray * arr, Int_t mode) const
2233 {
2234   //
2235   //sort tracks in array according mode criteria
2236   Int_t nseed = arr->GetEntriesFast();
2237   for (Int_t i=0; i<nseed; i++) {
2238     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2239     if (!pt) {
2240       continue;
2241     }
2242     pt->SetSort(mode);
2243   }
2244   arr->UnSort();
2245   arr->Sort();
2246 }
2247
2248
2249 void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t factor2, Int_t minimal)
2250 {
2251   //
2252   // Loop over all tracks and remove overlaped tracks (with lower quality)
2253   // Algorithm:
2254   //    1. Unsign clusters
2255   //    2. Sort tracks according quality
2256   //       Quality is defined by the number of cluster between first and last points
2257   //
2258   //    3. Loop over tracks - decreasing quality order
2259   //       a.) remove - If the fraction of shared cluster less than factor (1- n or 2)
2260   //       b.) remove - If the minimal number of clusters less than minimal and not ITS
2261   //       c.) if track accepted - sign clusters
2262   //
2263   //Called in - AliTPCtrackerMI::Clusters2Tracks()
2264   //          - AliTPCtrackerMI::PropagateBack()
2265   //          - AliTPCtrackerMI::RefitInward()
2266   //
2267
2268
2269   UnsignClusters();
2270   //
2271   Int_t nseed = arr->GetEntriesFast();
2272   Float_t * quality = new Float_t[nseed];
2273   Int_t   * indexes = new Int_t[nseed];
2274   Int_t good =0;
2275   //
2276   //
2277   for (Int_t i=0; i<nseed; i++) {
2278     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2279     if (!pt){
2280       quality[i]=-1;
2281       continue;
2282     }
2283     pt->UpdatePoints();    //select first last max dens points
2284     Float_t * points = pt->GetPoints();
2285     if (points[3]<0.8) quality[i] =-1;
2286     quality[i] = (points[2]-points[0])+pt->GetNumberOfClusters();
2287     //prefer high momenta tracks if overlaps
2288     quality[i] *= TMath::Sqrt(TMath::Abs(pt->Pt())+0.5);
2289   }
2290   TMath::Sort(nseed,quality,indexes);
2291   //
2292   //
2293   for (Int_t itrack=0; itrack<nseed; itrack++) {
2294     Int_t trackindex = indexes[itrack];
2295     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(trackindex);
2296     if (!pt) continue;
2297     //
2298     if (quality[trackindex]<0){
2299       if (pt) {
2300         delete arr->RemoveAt(trackindex);
2301       }
2302       else{
2303         arr->RemoveAt(trackindex);
2304       }
2305       continue;
2306     }
2307     //
2308     //
2309     Int_t first = Int_t(pt->GetPoints()[0]);
2310     Int_t last  = Int_t(pt->GetPoints()[2]);
2311     Double_t factor = (pt->GetBConstrain()) ? factor1: factor2;
2312     //
2313     Int_t found,foundable,shared;
2314     pt->GetClusterStatistic(first,last, found, foundable,shared,kFALSE); // better to get statistic in "high-dens"  region do't use full track as in line bellow
2315     //    pt->GetClusterStatistic(0,160, found, foundable,shared,kFALSE);
2316     Bool_t itsgold =kFALSE;
2317     if (pt->GetESD()){
2318       Int_t dummy[12];
2319       if (pt->GetESD()->GetITSclusters(dummy)>4) itsgold= kTRUE;
2320     }
2321     if (!itsgold){
2322       //
2323       if (Float_t(shared+1)/Float_t(found+1)>factor){
2324         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2325         delete arr->RemoveAt(trackindex);
2326         continue;
2327       }
2328       if (pt->GetNumberOfClusters()<50&&(found-0.5*shared)<minimal){  //remove short tracks
2329         if (pt->GetKinkIndexes()[0]!=0) continue;  //don't remove tracks  - part of the kinks
2330         delete arr->RemoveAt(trackindex);
2331         continue;
2332       }
2333     }
2334
2335     good++;
2336     //if (sharedfactor>0.4) continue;
2337     if (pt->GetKinkIndexes()[0]>0) continue;
2338     //Remove tracks with undefined properties - seems
2339     if (pt->GetSigmaY2()<kAlmost0) continue; // ? what is the origin ?
2340     //
2341     for (Int_t i=first; i<last; i++) {
2342       Int_t index=pt->GetClusterIndex2(i);
2343       // if (index<0 || index&0x8000 ) continue;
2344       if (index<0 || index&0x8000 ) continue;
2345       AliTPCclusterMI *c= pt->GetClusterPointer(i);
2346       if (!c) continue;
2347       c->Use(10);
2348     }
2349   }
2350   fNtracks = good;
2351   if (fDebug>0){
2352     Info("RemoveUsed","\n*****\nNumber of good tracks after shared removal\t%d\n",fNtracks);
2353   }
2354   delete []quality;
2355   delete []indexes;
2356 }
2357
2358 void AliTPCtrackerMI::UnsignClusters()
2359 {
2360   //
2361   // loop over all clusters and unsign them
2362   //
2363
2364   for (Int_t sec=0;sec<fkNIS;sec++){
2365     for (Int_t row=0;row<fInnerSec->GetNRows();row++){
2366       AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
2367       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++)
2368         //      if (cl[icl].IsUsed(10))
2369         cl[icl].Use(-1);
2370       cl = fInnerSec[sec][row].GetClusters2();
2371       for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++)
2372         //if (cl[icl].IsUsed(10))
2373           cl[icl].Use(-1);
2374     }
2375   }
2376
2377   for (Int_t sec=0;sec<fkNOS;sec++){
2378     for (Int_t row=0;row<fOuterSec->GetNRows();row++){
2379       AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
2380       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++)
2381         //if (cl[icl].IsUsed(10))
2382           cl[icl].Use(-1);
2383       cl = fOuterSec[sec][row].GetClusters2();
2384       for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++)
2385         //if (cl[icl].IsUsed(10))
2386         cl[icl].Use(-1);
2387     }
2388   }
2389
2390 }
2391
2392
2393
2394 void AliTPCtrackerMI::SignClusters(TObjArray * arr, Float_t fnumber, Float_t fdensity)
2395 {
2396   //
2397   //sign clusters to be "used"
2398   //
2399   // snumber and sdensity sign number of sigmas - bellow mean value to be accepted
2400   // loop over "primaries"
2401
2402   Float_t sumdens=0;
2403   Float_t sumdens2=0;
2404   Float_t sumn   =0;
2405   Float_t sumn2  =0;
2406   Float_t sumchi =0;
2407   Float_t sumchi2 =0;
2408
2409   Float_t sum    =0;
2410
2411   TStopwatch timer;
2412   timer.Start();
2413
2414   Int_t nseed = arr->GetEntriesFast();
2415   for (Int_t i=0; i<nseed; i++) {
2416     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2417     if (!pt) {
2418       continue;
2419     }
2420     if (!(pt->IsActive())) continue;
2421     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2422     if ( (dens>0.7) && (pt->GetNumberOfClusters()>70)){
2423       sumdens += dens;
2424       sumdens2+= dens*dens;
2425       sumn    += pt->GetNumberOfClusters();
2426       sumn2   += pt->GetNumberOfClusters()*pt->GetNumberOfClusters();
2427       Float_t chi2 = pt->GetChi2()/pt->GetNumberOfClusters();
2428       if (chi2>5) chi2=5;
2429       sumchi  +=chi2;
2430       sumchi2 +=chi2*chi2;
2431       sum++;
2432     }
2433   }
2434
2435   Float_t mdensity = 0.9;
2436   Float_t meann    = 130;
2437   Float_t meanchi  = 1;
2438   Float_t sdensity = 0.1;
2439   Float_t smeann    = 10;
2440   Float_t smeanchi  =0.4;
2441
2442
2443   if (sum>20){
2444     mdensity = sumdens/sum;
2445     meann    = sumn/sum;
2446     meanchi  = sumchi/sum;
2447     //
2448     sdensity = sumdens2/sum-mdensity*mdensity;
2449     if (sdensity >= 0)
2450        sdensity = TMath::Sqrt(sdensity);
2451     else
2452        sdensity = 0.1;
2453     //
2454     smeann   = sumn2/sum-meann*meann;
2455     if (smeann >= 0)
2456       smeann   = TMath::Sqrt(smeann);
2457     else
2458       smeann   = 10;
2459     //
2460     smeanchi = sumchi2/sum - meanchi*meanchi;
2461     if (smeanchi >= 0)
2462       smeanchi = TMath::Sqrt(smeanchi);
2463     else
2464       smeanchi = 0.4;
2465   }
2466
2467
2468   //REMOVE  SHORT DELTAS or tracks going out of sensitive volume of TPC
2469   //
2470   for (Int_t i=0; i<nseed; i++) {
2471     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2472     if (!pt) {
2473       continue;
2474     }
2475     if (pt->GetBSigned()) continue;
2476     if (pt->GetBConstrain()) continue;
2477     //if (!(pt->IsActive())) continue;
2478     /*
2479     Int_t found,foundable,shared;
2480     pt->GetClusterStatistic(0,160,found, foundable,shared);
2481     if (shared/float(found)>0.3) {
2482       if (shared/float(found)>0.9 ){
2483         //delete arr->RemoveAt(i);
2484       }
2485       continue;
2486     }
2487     */
2488     Bool_t isok =kFALSE;
2489     if ( (pt->GetNShared()/pt->GetNumberOfClusters()<0.5) &&pt->GetNumberOfClusters()>60)
2490       isok = kTRUE;
2491     if ((TMath::Abs(1/pt->GetC())<100.) && (pt->GetNShared()/pt->GetNumberOfClusters()<0.7))
2492       isok =kTRUE;
2493     if  (TMath::Abs(pt->GetZ()/pt->GetX())>1.1)
2494       isok =kTRUE;
2495     if ( (TMath::Abs(pt->GetSnp()>0.7) && pt->GetD(0,0)>60.))
2496       isok =kTRUE;
2497
2498     if (isok)
2499       for (Int_t i=0; i<160; i++) {
2500         Int_t index=pt->GetClusterIndex2(i);
2501         if (index<0) continue;
2502         AliTPCclusterMI *c= pt->GetClusterPointer(i);
2503         if (!c) continue;
2504         //if (!(c->IsUsed(10))) c->Use();
2505         c->Use(10);
2506       }
2507   }
2508
2509
2510   //
2511   Double_t maxchi  = meanchi+2.*smeanchi;
2512
2513   for (Int_t i=0; i<nseed; i++) {
2514     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2515     if (!pt) {
2516       continue;
2517     }
2518     //if (!(pt->IsActive())) continue;
2519     if (pt->GetBSigned()) continue;
2520     Double_t chi     = pt->GetChi2()/pt->GetNumberOfClusters();
2521     if (chi>maxchi) continue;
2522
2523     Float_t bfactor=1;
2524     Float_t dens = pt->GetNumberOfClusters()/Float_t(pt->GetNFoundable());
2525
2526     //sign only tracks with enoug big density at the beginning
2527
2528     if ((pt->GetDensityFirst(40)<0.75) && pt->GetNumberOfClusters()<meann) continue;
2529
2530
2531     Double_t mindens = TMath::Max(double(mdensity-sdensity*fdensity*bfactor),0.65);
2532     Double_t minn    = TMath::Max(Int_t(meann-fnumber*smeann*bfactor),50);
2533
2534     //    if (pt->fBConstrain) mindens = TMath::Max(mdensity-sdensity*fdensity*bfactor,0.65);
2535     if ( (pt->GetRemoval()==10) && (pt->GetSnp()>0.8)&&(dens>mindens))
2536       minn=0;
2537
2538     if ((dens>mindens && pt->GetNumberOfClusters()>minn) && chi<maxchi ){
2539       //Int_t noc=pt->GetNumberOfClusters();
2540       pt->SetBSigned(kTRUE);
2541       for (Int_t i=0; i<160; i++) {
2542
2543         Int_t index=pt->GetClusterIndex2(i);
2544         if (index<0) continue;
2545         AliTPCclusterMI *c= pt->GetClusterPointer(i);
2546         if (!c) continue;
2547         //      if (!(c->IsUsed(10))) c->Use();
2548         c->Use(10);
2549       }
2550     }
2551   }
2552   //  gLastCheck = nseed;
2553   //  arr->Compress();
2554   if (fDebug>0){
2555     timer.Print();
2556   }
2557 }
2558
2559
2560 void  AliTPCtrackerMI::StopNotActive(TObjArray * arr, Int_t row0, Float_t th0, Float_t th1, Float_t th2) const
2561 {
2562   // stop not active tracks
2563   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2564   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2565   Int_t nseed = arr->GetEntriesFast();
2566   //
2567   for (Int_t i=0; i<nseed; i++) {
2568     AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
2569     if (!pt) {
2570       continue;
2571     }
2572     if (!(pt->IsActive())) continue;
2573     StopNotActive(pt,row0,th0, th1,th2);
2574   }
2575 }
2576
2577
2578
2579 void  AliTPCtrackerMI::StopNotActive(AliTPCseed * seed, Int_t row0, Float_t th0, Float_t th1,
2580  Float_t th2) const
2581 {
2582   // stop not active tracks
2583   // take th1 as threshold for number of founded to number of foundable on last 10 active rows
2584   // take th2 as threshold for number of founded to number of foundable on last 20 active rows
2585   Int_t sumgood1  = 0;
2586   Int_t sumgood2  = 0;
2587   Int_t foundable = 0;
2588   Int_t maxindex = seed->GetLastPoint();  //last foundable row
2589   if (seed->GetNFoundable()*th0 > seed->GetNumberOfClusters()) {
2590     seed->Desactivate(10) ;
2591     return;
2592   }
2593
2594   for (Int_t i=row0; i<maxindex; i++){
2595     Int_t index = seed->GetClusterIndex2(i);
2596     if (index!=-1) foundable++;
2597     //if (!c) continue;
2598     if (foundable<=30) sumgood1++;
2599     if (foundable<=50) {
2600       sumgood2++;
2601     }
2602     else{
2603       break;
2604     }
2605   }
2606   if (foundable>=30.){
2607      if (sumgood1<(th1*30.)) seed->Desactivate(10);
2608   }
2609   if (foundable>=50)
2610     if (sumgood2<(th2*50.)) seed->Desactivate(10);
2611 }
2612
2613
2614 Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
2615 {
2616   //
2617   // back propagation of ESD tracks
2618   //
2619   //return 0;
2620   const Int_t kMaxFriendTracks=2000;
2621   fEvent = event;
2622   ReadSeeds(event,2);
2623   fIteration=2;
2624   //PrepareForProlongation(fSeeds,1);
2625   PropagateForward2(fSeeds);
2626   RemoveUsed2(fSeeds,0.4,0.4,20);
2627
2628   TObjArray arraySeed(fSeeds->GetEntries());
2629   for (Int_t i=0;i<fSeeds->GetEntries();i++) {
2630     arraySeed.AddAt(fSeeds->At(i),i);
2631   }
2632   SignShared(&arraySeed);
2633   //  FindCurling(fSeeds, event,2); // find multi found tracks
2634   FindSplitted(fSeeds, event,2); // find multi found tracks
2635
2636   Int_t ntracks=0;
2637   Int_t nseed = fSeeds->GetEntriesFast();
2638   for (Int_t i=0;i<nseed;i++){
2639     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2640     if (!seed) continue;
2641     if (seed->GetKinkIndex(0)>0)  UpdateKinkQualityD(seed);  // update quality informations for kinks
2642
2643     seed->PropagateTo(fParam->GetInnerRadiusLow());
2644     seed->UpdatePoints();
2645     AddCovariance(seed);
2646     MakeBitmaps(seed);
2647     AliESDtrack *esd=event->GetTrack(i);
2648     seed->CookdEdx(0.02,0.6);
2649     CookLabel(seed,0.1); //For comparison only
2650     esd->SetTPCClusterMap(seed->GetClusterMap());
2651     esd->SetTPCSharedMap(seed->GetSharedMap());
2652     //
2653     if (AliTPCReconstructor::StreamLevel()>1 && seed!=0&&esd!=0) {
2654       TTreeSRedirector &cstream = *fDebugStreamer;
2655       cstream<<"Crefit"<<
2656         "Esd.="<<esd<<
2657         "Track.="<<seed<<
2658         "\n";
2659     }
2660
2661     if (seed->GetNumberOfClusters()>15){
2662       esd->UpdateTrackParams(seed,AliESDtrack::kTPCrefit);
2663       esd->SetTPCPoints(seed->GetPoints());
2664       esd->SetTPCPointsF(seed->GetNFoundable());
2665       Int_t ndedx   = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2666       Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2667       Float_t dedx  = seed->GetdEdx();
2668       esd->SetTPCsignal(dedx, sdedx, ndedx);
2669       //
2670       // add seed to the esd track in Calib level
2671       //
2672       Bool_t storeFriend = gRandom->Rndm()<(kMaxFriendTracks)/Float_t(nseed);
2673       if (AliTPCReconstructor::StreamLevel()>0 &&storeFriend){
2674         AliTPCseed * seedCopy = new AliTPCseed(*seed, kTRUE);
2675         esd->AddCalibObject(seedCopy);
2676       }
2677       ntracks++;
2678     }
2679     else{
2680       //printf("problem\n");
2681     }
2682   }
2683   //FindKinks(fSeeds,event);
2684   Info("RefitInward","Number of refitted tracks %d",ntracks);
2685   return 0;
2686 }
2687
2688
2689 Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
2690 {
2691   //
2692   // back propagation of ESD tracks
2693   //
2694
2695   fEvent = event;
2696   fIteration = 1;
2697   ReadSeeds(event,1);
2698   PropagateBack(fSeeds);
2699   RemoveUsed2(fSeeds,0.4,0.4,20);
2700   //FindCurling(fSeeds, fEvent,1);
2701   FindSplitted(fSeeds, event,1); // find multi found tracks
2702
2703   //
2704   Int_t nseed = fSeeds->GetEntriesFast();
2705   Int_t ntracks=0;
2706   for (Int_t i=0;i<nseed;i++){
2707     AliTPCseed * seed = (AliTPCseed*) fSeeds->UncheckedAt(i);
2708     if (!seed) continue;
2709     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
2710     seed->UpdatePoints();
2711     AddCovariance(seed);
2712     AliESDtrack *esd=event->GetTrack(i);
2713     seed->CookdEdx(0.02,0.6);
2714     CookLabel(seed,0.1); //For comparison only
2715     if (seed->GetNumberOfClusters()>15){
2716       esd->UpdateTrackParams(seed,AliESDtrack::kTPCout);
2717       esd->SetTPCPoints(seed->GetPoints());
2718       esd->SetTPCPointsF(seed->GetNFoundable());
2719       Int_t ndedx   = seed->GetNCDEDX(0)+seed->GetNCDEDX(1)+seed->GetNCDEDX(2)+seed->GetNCDEDX(3);
2720       Float_t sdedx = (seed->GetSDEDX(0)+seed->GetSDEDX(1)+seed->GetSDEDX(2)+seed->GetSDEDX(3))*0.25;
2721       Float_t dedx  = seed->GetdEdx();
2722       esd->SetTPCsignal(dedx, sdedx, ndedx);
2723       ntracks++;
2724       Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
2725       // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number
2726       if (AliTPCReconstructor::StreamLevel()>1 && esd) {
2727         (*fDebugStreamer)<<"Cback"<<
2728           "Tr0.="<<seed<<
2729           "esd.="<<esd<<
2730           "EventNrInFile="<<eventnumber<<
2731           "\n";
2732       }
2733     }
2734   }
2735   //FindKinks(fSeeds,event);
2736   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
2737   fEvent =0;
2738
2739   return 0;
2740 }
2741
2742
2743 void AliTPCtrackerMI::DeleteSeeds()
2744 {
2745   //
2746   //delete Seeds
2747
2748   Int_t nseed = fSeeds->GetEntriesFast();
2749   for (Int_t i=0;i<nseed;i++){
2750     AliTPCseed * seed = (AliTPCseed*)fSeeds->At(i);
2751     if (seed) delete fSeeds->RemoveAt(i);
2752   }
2753   delete fSeeds;
2754
2755   fSeeds =0;
2756 }
2757
2758 void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
2759 {
2760   //
2761   //read seeds from the event
2762
2763   Int_t nentr=event->GetNumberOfTracks();
2764   if (fDebug>0){
2765     Info("ReadSeeds", "Number of ESD tracks: %d\n", nentr);
2766   }
2767   if (fSeeds)
2768     DeleteSeeds();
2769   if (!fSeeds){
2770     fSeeds = new TObjArray(nentr);
2771   }
2772   UnsignClusters();
2773   //  Int_t ntrk=0;
2774   for (Int_t i=0; i<nentr; i++) {
2775     AliESDtrack *esd=event->GetTrack(i);
2776     ULong_t status=esd->GetStatus();
2777     if (!(status&AliESDtrack::kTPCin)) continue;
2778     AliTPCtrack t(*esd);
2779     t.SetNumberOfClusters(0);
2780     //    AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
2781     AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
2782     seed->SetUniqueID(esd->GetID());
2783     AddCovariance(seed);   //add systematic ucertainty
2784     for (Int_t ikink=0;ikink<3;ikink++) {
2785       Int_t index = esd->GetKinkIndex(ikink);
2786       seed->GetKinkIndexes()[ikink] = index;
2787       if (index==0) continue;
2788       index = TMath::Abs(index);
2789       AliESDkink * kink = fEvent->GetKink(index-1);
2790       if (kink&&esd->GetKinkIndex(ikink)<0){
2791         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,2);
2792         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,0);
2793       }
2794       if (kink&&esd->GetKinkIndex(ikink)>0){
2795         if ((status & AliESDtrack::kTRDrefit) != 0) kink->SetStatus(1,6);
2796         if ((status & AliESDtrack::kITSout) != 0)   kink->SetStatus(1,4);
2797       }
2798
2799     }
2800     if (((status&AliESDtrack::kITSout)==0)&&(direction==1)) seed->ResetCovariance(10.);
2801     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance(10.);
2802     if ( direction ==2 && ((status & AliESDtrack::kTPCout) == 0) ) {
2803       fSeeds->AddAt(0,i);
2804       delete seed;
2805       continue;
2806     }
2807     if ( direction ==2 &&(status & AliESDtrack::kTRDrefit) > 0 )  {
2808       Double_t par0[5],par1[5],alpha,x;
2809       esd->GetInnerExternalParameters(alpha,x,par0);
2810       esd->GetExternalParameters(x,par1);
2811       Double_t delta1 = TMath::Abs(par0[4]-par1[4])/(0.000000001+TMath::Abs(par0[4]+par1[4]));
2812       Double_t delta2 = TMath::Abs(par0[3]-par1[3]);
2813       Double_t trdchi2=0;
2814       if (esd->GetTRDncls()>0) trdchi2 = esd->GetTRDchi2()/esd->GetTRDncls();
2815       //reset covariance if suspicious
2816       if ( (delta1>0.1) || (delta2>0.006) ||trdchi2>7.)
2817         seed->ResetCovariance(10.);
2818     }
2819
2820     //
2821     //
2822     // rotate to the local coordinate system
2823     //
2824     fSectors=fInnerSec; fN=fkNIS;
2825     Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
2826     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
2827     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
2828     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
2829     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
2830     if (alpha<-TMath::Pi()) alpha += 2*TMath::Pi();
2831     if (alpha>=TMath::Pi()) alpha -= 2*TMath::Pi();
2832     alpha-=seed->GetAlpha();
2833     if (!seed->Rotate(alpha)) {
2834       delete seed;
2835       continue;
2836     }
2837     seed->SetESD(esd);
2838     // sign clusters
2839     if (esd->GetKinkIndex(0)<=0){
2840       for (Int_t irow=0;irow<160;irow++){
2841         Int_t index = seed->GetClusterIndex2(irow);
2842         if (index>0){
2843           //
2844           AliTPCclusterMI * cl = GetClusterMI(index);
2845           seed->SetClusterPointer(irow,cl);
2846           if (cl){
2847             if ((index & 0x8000)==0){
2848               cl->Use(10);  // accepted cluster
2849             }else{
2850               cl->Use(6);   // close cluster not accepted
2851             }
2852           }else{
2853             Info("ReadSeeds","Not found cluster");
2854           }
2855         }
2856       }
2857     }
2858     fSeeds->AddAt(seed,i);
2859   }
2860 }
2861
2862
2863
2864 //_____________________________________________________________________________
2865 void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
2866                                  Float_t deltay, Int_t ddsec) {
2867   //-----------------------------------------------------------------
2868   // This function creates track seeds.
2869   // SEEDING WITH VERTEX CONSTRAIN
2870   //-----------------------------------------------------------------
2871   // cuts[0]   - fP4 cut
2872   // cuts[1]   - tan(phi)  cut
2873   // cuts[2]   - zvertex cut
2874   // cuts[3]   - fP3 cut
2875   Int_t nin0  = 0;
2876   Int_t nin1  = 0;
2877   Int_t nin2  = 0;
2878   Int_t nin   = 0;
2879   Int_t nout1 = 0;
2880   Int_t nout2 = 0;
2881
2882   Double_t x[5], c[15];
2883   //  Int_t di = i1-i2;
2884   //
2885   AliTPCseed * seed = new AliTPCseed();
2886   Double_t alpha=fSectors->GetAlpha(), shift=fSectors->GetAlphaShift();
2887   Double_t cs=cos(alpha), sn=sin(alpha);
2888   //
2889   //  Double_t x1 =fOuterSec->GetX(i1);
2890   //Double_t xx2=fOuterSec->GetX(i2);
2891
2892   Double_t x1 =GetXrow(i1);
2893   Double_t xx2=GetXrow(i2);
2894
2895   Double_t x3=GetX(), y3=GetY(), z3=GetZ();
2896
2897   Int_t imiddle = (i2+i1)/2;    //middle pad row index
2898   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
2899   const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
2900   //
2901   Int_t ns =sec;
2902
2903   const AliTPCtrackerRow& kr1=GetRow(ns,i1);
2904   Double_t ymax  = GetMaxY(i1)-kr1.GetDeadZone()-1.5;
2905   Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;
2906
2907   //
2908   // change cut on curvature if it can't reach this layer
2909   // maximal curvature set to reach it
2910   Double_t dvertexmax  = TMath::Sqrt((x1-x3)*(x1-x3)+(ymax+5-y3)*(ymax+5-y3));
2911   if (dvertexmax*0.5*cuts[0]>0.85){
2912     cuts[0] = 0.85/(dvertexmax*0.5+1.);
2913   }
2914   Double_t r2min = 1/(cuts[0]*cuts[0]);  //minimal square of radius given by cut
2915
2916   //  Int_t ddsec = 1;
2917   if (deltay>0) ddsec = 0;
2918   // loop over clusters
2919   for (Int_t is=0; is < kr1; is++) {
2920     //
2921     if (kr1[is]->IsUsed(10)) continue;
2922     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
2923     //if (TMath::Abs(y1)>ymax) continue;
2924
2925     if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
2926
2927     // find possible directions
2928     Float_t anglez = (z1-z3)/(x1-x3);
2929     Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
2930     //
2931     //
2932     //find   rotation angles relative to line given by vertex and point 1
2933     Double_t dvertex2 = (x1-x3)*(x1-x3)+(y1-y3)*(y1-y3);
2934     Double_t dvertex  = TMath::Sqrt(dvertex2);
2935     Double_t angle13  = TMath::ATan((y1-y3)/(x1-x3));
2936     Double_t cs13     = cos(-angle13), sn13 = sin(-angle13);
2937
2938     //
2939     // loop over 2 sectors
2940     Int_t dsec1=-ddsec;
2941     Int_t dsec2= ddsec;
2942     if (y1<0)  dsec2= 0;
2943     if (y1>0)  dsec1= 0;
2944
2945     Double_t dddz1=0;  // direction of delta inclination in z axis
2946     Double_t dddz2=0;
2947     if ( (z1-z3)>0)
2948       dddz1 =1;
2949     else
2950       dddz2 =1;
2951     //
2952     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
2953       Int_t sec2 = sec + dsec;
2954       //
2955       //      AliTPCtrackerRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
2956       //AliTPCtrackerRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
2957       AliTPCtrackerRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
2958       AliTPCtrackerRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
2959       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
2960       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
2961
2962       // rotation angles to p1-p3
2963       Double_t cs13r     = cos(-angle13+dsec*alpha)/dvertex, sn13r = sin(-angle13+dsec*alpha)/dvertex;
2964       Double_t x2,   y2,   z2;
2965       //
2966       //      Double_t dymax = maxangle*TMath::Abs(x1-xx2);
2967
2968       //
2969       Double_t dxx0 =  (xx2-x3)*cs13r;
2970       Double_t dyy0 =  (xx2-x3)*sn13r;
2971       for (Int_t js=index1; js < index2; js++) {
2972         const AliTPCclusterMI *kcl = kr2[js];
2973         if (kcl->IsUsed(10)) continue;
2974         //
2975         //calcutate parameters
2976         //
2977         Double_t yy0 =  dyy0 +(kcl->GetY()-y3)*cs13r;
2978         // stright track
2979         if (TMath::Abs(yy0)<0.000001) continue;
2980         Double_t xx0 =  dxx0 -(kcl->GetY()-y3)*sn13r;
2981         Double_t y0  =  0.5*(xx0*xx0+yy0*yy0-xx0)/yy0;
2982         Double_t r02 = (0.25+y0*y0)*dvertex2;
2983         //curvature (radius) cut
2984         if (r02<r2min) continue;
2985
2986         nin0++;
2987         //
2988         Double_t c0  = 1/TMath::Sqrt(r02);
2989         if (yy0>0) c0*=-1.;
2990
2991
2992         //Double_t dfi0   = 2.*TMath::ASin(dvertex*c0*0.5);
2993         //Double_t dfi1   = 2.*TMath::ASin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2994         Double_t dfi0   = 2.*AliTPCFastMath::FastAsin(dvertex*c0*0.5);
2995         Double_t dfi1   = 2.*AliTPCFastMath::FastAsin(TMath::Sqrt(yy0*yy0+(1-xx0)*(1-xx0))*dvertex*c0*0.5);
2996         //
2997         //
2998         Double_t z0  =  kcl->GetZ();
2999         Double_t zzzz2    = z1-(z1-z3)*dfi1/dfi0;
3000         if (TMath::Abs(zzzz2-z0)>0.5) continue;
3001         nin1++;
3002         //
3003         Double_t dip    = (z1-z0)*c0/dfi1;
3004         Double_t x0 = (0.5*cs13+y0*sn13)*dvertex*c0;
3005         //
3006         y2 = kcl->GetY();
3007         if (dsec==0){
3008           x2 = xx2;
3009           z2 = kcl->GetZ();
3010         }
3011         else
3012           {
3013             // rotation
3014             z2 = kcl->GetZ();
3015             x2= xx2*cs-y2*sn*dsec;
3016             y2=+xx2*sn*dsec+y2*cs;
3017           }
3018
3019         x[0] = y1;
3020         x[1] = z1;
3021         x[2] = x0;
3022         x[3] = dip;
3023         x[4] = c0;
3024         //
3025         //
3026         // do we have cluster at the middle ?
3027         Double_t ym,zm;
3028         GetProlongation(x1,xm,x,ym,zm);
3029         UInt_t dummy;
3030         AliTPCclusterMI * cm=0;
3031         if (TMath::Abs(ym)-ymaxm<0){
3032           cm = krm.FindNearest2(ym,zm,1.0,0.6,dummy);
3033           if ((!cm) || (cm->IsUsed(10))) {
3034             continue;
3035           }
3036         }
3037         else{
3038           // rotate y1 to system 0
3039           // get state vector in rotated system
3040           Double_t yr1  = (-0.5*sn13+y0*cs13)*dvertex*c0;
3041           Double_t xr2  =  x0*cs+yr1*sn*dsec;
3042           Double_t xr[5]={kcl->GetY(),kcl->GetZ(), xr2, dip, c0};
3043           //
3044           GetProlongation(xx2,xm,xr,ym,zm);
3045           if (TMath::Abs(ym)-ymaxm<0){
3046             cm = kr2m.FindNearest2(ym,zm,1.0,0.6,dummy);
3047             if ((!cm) || (cm->IsUsed(10))) {
3048               continue;
3049             }
3050           }
3051         }
3052
3053
3054         Double_t dym = 0;
3055         Double_t dzm = 0;
3056         if (cm){
3057           dym = ym - cm->GetY();
3058           dzm = zm - cm->GetZ();
3059         }
3060         nin2++;
3061
3062
3063         //
3064         //
3065         Double_t sy1=kr1[is]->GetSigmaY2()*2., sz1=kr1[is]->GetSigmaZ2()*2.;
3066         Double_t sy2=kcl->GetSigmaY2()*2.,     sz2=kcl->GetSigmaZ2()*2.;
3067         //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
3068         Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
3069         //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
3070
3071         Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3072         Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3073         Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3074         Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3075         Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3076         Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3077
3078         Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3079         Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3080         Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3081         Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3082
3083         c[0]=sy1;
3084         c[1]=0.;       c[2]=sz1;
3085         c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3086         c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3087                        c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3088         c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3089         c[13]=f30*sy1*f40+f32*sy2*f42;
3090         c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3091
3092         //      if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3093
3094         UInt_t index=kr1.GetIndex(is);
3095         seed->~AliTPCseed(); // this does not set the pointer to 0...
3096         AliTPCseed *track=new(seed) AliTPCseed(x1, ns*alpha+shift, x, c, index);
3097
3098         track->SetIsSeeding(kTRUE);
3099         track->SetSeed1(i1);
3100         track->SetSeed2(i2);
3101         track->SetSeedType(3);
3102
3103
3104         //if (dsec==0) {
3105           FollowProlongation(*track, (i1+i2)/2,1);
3106           Int_t foundable,found,shared;
3107           track->GetClusterStatistic((i1+i2)/2,i1, found, foundable, shared, kTRUE);
3108           if ((found<0.55*foundable)  || shared>0.5*found || (track->GetSigmaY2()+track->GetSigmaZ2())>0.5){
3109             seed->Reset();
3110             seed->~AliTPCseed();
3111             continue;
3112           }
3113           //}
3114
3115         nin++;
3116         FollowProlongation(*track, i2,1);
3117
3118
3119         //Int_t rc = 1;
3120         track->SetBConstrain(1);
3121         //      track->fLastPoint = i1+fInnerSec->GetNRows();  // first cluster in track position
3122         track->SetLastPoint(i1);  // first cluster in track position
3123         track->SetFirstPoint(track->GetLastPoint());
3124
3125         if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3126             track->GetNumberOfClusters() < track->GetNFoundable()*0.6 ||
3127             track->GetNShared()>0.4*track->GetNumberOfClusters() ) {
3128           seed->Reset();
3129           seed->~AliTPCseed();
3130           continue;
3131         }
3132         nout1++;
3133         // Z VERTEX CONDITION
3134         Double_t zv, bz=GetBz();
3135         if ( !track->GetZAt(0.,bz,zv) ) continue;
3136         if (TMath::Abs(zv-z3)>cuts[2]) {
3137           FollowProlongation(*track, TMath::Max(i2-20,0));
3138           if ( !track->GetZAt(0.,bz,zv) ) continue;
3139           if (TMath::Abs(zv-z3)>cuts[2]){
3140             FollowProlongation(*track, TMath::Max(i2-40,0));
3141             if ( !track->GetZAt(0.,bz,zv) ) continue;
3142             if (TMath::Abs(zv-z3)>cuts[2] &&(track->GetNumberOfClusters() > track->GetNFoundable()*0.7)){
3143               // make seed without constrain
3144               AliTPCseed * track2 = MakeSeed(track,0.2,0.5,1.);
3145               FollowProlongation(*track2, i2,1);
3146               track2->SetBConstrain(kFALSE);
3147               track2->SetSeedType(1);
3148               arr->AddLast(track2);
3149               seed->Reset();
3150               seed->~AliTPCseed();
3151               continue;
3152             }
3153             else{
3154               seed->Reset();
3155               seed->~AliTPCseed();
3156               continue;
3157
3158             }
3159           }
3160         }
3161
3162         track->SetSeedType(0);
3163         arr->AddLast(track);
3164         seed = new AliTPCseed;
3165         nout2++;
3166         // don't consider other combinations
3167         if (track->GetNumberOfClusters() > track->GetNFoundable()*0.8)
3168           break;
3169       }
3170     }
3171   }
3172   if (fDebug>3){
3173     Info("MakeSeeds3","\nSeeding statistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2);
3174   }
3175   delete seed;
3176 }
3177
3178
3179 void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,  Float_t cuts[4],
3180                                  Float_t deltay) {
3181
3182
3183
3184   //-----------------------------------------------------------------
3185   // This function creates track seeds.
3186   //-----------------------------------------------------------------
3187   // cuts[0]   - fP4 cut
3188   // cuts[1]   - tan(phi)  cut
3189   // cuts[2]   - zvertex cut
3190   // cuts[3]   - fP3 cut
3191
3192
3193   Int_t nin0  = 0;
3194   Int_t nin1  = 0;
3195   Int_t nin2  = 0;
3196   Int_t nin   = 0;
3197   Int_t nout1 = 0;
3198   Int_t nout2 = 0;
3199   Int_t nout3 =0;
3200   Double_t x[5], c[15];
3201   //
3202   // make temporary seed
3203   AliTPCseed * seed = new AliTPCseed;
3204   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3205   //  Double_t cs=cos(alpha), sn=sin(alpha);
3206   //
3207   //
3208
3209   // first 3 padrows
3210   Double_t x1 = GetXrow(i1-1);
3211   const    AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
3212   Double_t y1max  = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;
3213   //
3214   Double_t x1p = GetXrow(i1);
3215   const    AliTPCtrackerRow& kr1p=GetRow(sec,i1);
3216   //
3217   Double_t x1m = GetXrow(i1-2);
3218   const    AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
3219
3220   //
3221   //last 3 padrow for seeding
3222   AliTPCtrackerRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
3223   Double_t    x3   =  GetXrow(i1-7);
3224   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;
3225   //
3226   AliTPCtrackerRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
3227   Double_t    x3p   = GetXrow(i1-6);
3228   //
3229   AliTPCtrackerRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
3230   Double_t    x3m   = GetXrow(i1-8);
3231
3232   //
3233   //
3234   // middle padrow
3235   Int_t im = i1-4;                           //middle pad row index
3236   Double_t xm         = GetXrow(im);         // radius of middle pad-row
3237   const AliTPCtrackerRow& krm=GetRow(sec,im);   //middle pad -row
3238   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;
3239   //
3240   //
3241   Double_t deltax  = x1-x3;
3242   Double_t dymax   = deltax*cuts[1];
3243   Double_t dzmax   = deltax*cuts[3];
3244   //
3245   // loop over clusters
3246   for (Int_t is=0; is < kr1; is++) {
3247     //
3248     if (kr1[is]->IsUsed(10)) continue;
3249     Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
3250     //
3251     if (deltay>0 && TMath::Abs(y1max-TMath::Abs(y1))> deltay ) continue;  // seed only at the edge
3252     //
3253     Int_t  index1 = TMath::Max(kr3.Find(z1-dzmax)-1,0);
3254     Int_t  index2 = TMath::Min(kr3.Find(z1+dzmax)+1,kr3);
3255     //
3256     Double_t y3,   z3;
3257     //
3258     //
3259     UInt_t index;
3260     for (Int_t js=index1; js < index2; js++) {
3261       const AliTPCclusterMI *kcl = kr3[js];
3262       if (kcl->IsUsed(10)) continue;
3263       y3 = kcl->GetY();
3264       // apply angular cuts
3265       if (TMath::Abs(y1-y3)>dymax) continue;
3266       x3 = x3;
3267       z3 = kcl->GetZ();
3268       if (TMath::Abs(z1-z3)>dzmax) continue;
3269       //
3270       Double_t angley = (y1-y3)/(x1-x3);
3271       Double_t anglez = (z1-z3)/(x1-x3);
3272       //
3273       Double_t erry = TMath::Abs(angley)*(x1-x1m)*0.5+0.5;
3274       Double_t errz = TMath::Abs(anglez)*(x1-x1m)*0.5+0.5;
3275       //
3276       Double_t yyym = angley*(xm-x1)+y1;
3277       Double_t zzzm = anglez*(xm-x1)+z1;
3278
3279       const AliTPCclusterMI *kcm = krm.FindNearest2(yyym,zzzm,erry,errz,index);
3280       if (!kcm) continue;
3281       if (kcm->IsUsed(10)) continue;
3282
3283       erry = TMath::Abs(angley)*(x1-x1m)*0.4+0.5;
3284       errz = TMath::Abs(anglez)*(x1-x1m)*0.4+0.5;
3285       //
3286       //
3287       //
3288       Int_t used  =0;
3289       Int_t found =0;
3290       //
3291       // look around first
3292       const AliTPCclusterMI *kc1m = kr1m.FindNearest2(angley*(x1m-x1)+y1,
3293                                                       anglez*(x1m-x1)+z1,
3294                                                       erry,errz,index);
3295       //
3296       if (kc1m){
3297         found++;
3298         if (kc1m->IsUsed(10)) used++;
3299       }
3300       const AliTPCclusterMI *kc1p = kr1p.FindNearest2(angley*(x1p-x1)+y1,
3301                                                       anglez*(x1p-x1)+z1,
3302                                                       erry,errz,index);
3303       //
3304       if (kc1p){
3305         found++;
3306         if (kc1p->IsUsed(10)) used++;
3307       }
3308       if (used>1)  continue;
3309       if (found<1) continue;
3310
3311       //
3312       // look around last
3313       const AliTPCclusterMI *kc3m = kr3m.FindNearest2(angley*(x3m-x3)+y3,
3314                                                       anglez*(x3m-x3)+z3,
3315                                                       erry,errz,index);
3316       //
3317       if (kc3m){
3318         found++;
3319         if (kc3m->IsUsed(10)) used++;
3320       }
3321       else
3322         continue;
3323       const AliTPCclusterMI *kc3p = kr3p.FindNearest2(angley*(x3p-x3)+y3,
3324                                                       anglez*(x3p-x3)+z3,
3325                                                       erry,errz,index);
3326       //
3327       if (kc3p){
3328         found++;
3329         if (kc3p->IsUsed(10)) used++;
3330       }
3331       else
3332         continue;
3333       if (used>1)  continue;
3334       if (found<3) continue;
3335       //
3336       Double_t x2,y2,z2;
3337       x2 = xm;
3338       y2 = kcm->GetY();
3339       z2 = kcm->GetZ();
3340       //
3341
3342       x[0]=y1;
3343       x[1]=z1;
3344       x[4]=F1(x1,y1,x2,y2,x3,y3);
3345       //if (TMath::Abs(x[4]) >= cuts[0]) continue;
3346       nin0++;
3347       //
3348       x[2]=F2(x1,y1,x2,y2,x3,y3);
3349       nin1++;
3350       //
3351       x[3]=F3n(x1,y1,x2,y2,z1,z2,x[4]);
3352       //if (TMath::Abs(x[3]) > cuts[3]) continue;
3353       nin2++;
3354       //
3355       //
3356       Double_t sy1=0.1,  sz1=0.1;
3357       Double_t sy2=0.1,  sz2=0.1;
3358       Double_t sy3=0.1,  sy=0.1, sz=0.1;
3359
3360       Double_t f40=(F1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
3361       Double_t f42=(F1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
3362       Double_t f43=(F1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
3363       Double_t f20=(F2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
3364       Double_t f22=(F2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
3365       Double_t f23=(F2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
3366
3367       Double_t f30=(F3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
3368       Double_t f31=(F3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
3369       Double_t f32=(F3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
3370       Double_t f34=(F3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
3371
3372       c[0]=sy1;
3373       c[1]=0.;       c[2]=sz1;
3374       c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
3375       c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
3376       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
3377       c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
3378       c[13]=f30*sy1*f40+f32*sy2*f42;
3379       c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
3380
3381       //        if (!BuildSeed(kr1[is],kcl,0,x1,x2,x3,x,c)) continue;
3382
3383       UInt_t index=kr1.GetIndex(is);
3384       seed->~AliTPCseed();
3385       AliTPCseed *track=new(seed) AliTPCseed(x1, sec*alpha+shift, x, c, index);
3386
3387       track->SetIsSeeding(kTRUE);
3388
3389       nin++;
3390       FollowProlongation(*track, i1-7,1);
3391       if (track->GetNumberOfClusters() < track->GetNFoundable()*0.75 ||
3392           track->GetNShared()>0.6*track->GetNumberOfClusters() || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.6){
3393         seed->Reset();
3394         seed->~AliTPCseed();
3395         continue;
3396       }
3397       nout1++;
3398       nout2++;
3399       //Int_t rc = 1;
3400       FollowProlongation(*track, i2,1);
3401       track->SetBConstrain(0);
3402       track->SetLastPoint(i1+fInnerSec->GetNRows());  // first cluster in track position
3403       track->SetFirstPoint(track->GetLastPoint());
3404
3405       if (track->GetNumberOfClusters()<(i1-i2)*0.5 ||
3406           track->GetNumberOfClusters()<track->GetNFoundable()*0.7 ||
3407           track->GetNShared()>2. || track->GetChi2()/track->GetNumberOfClusters()>6 || ( track->GetSigmaY2()+ track->GetSigmaZ2())>0.5 ) {
3408         seed->Reset();
3409         seed->~AliTPCseed();
3410         continue;
3411       }
3412
3413       {
3414         FollowProlongation(*track, TMath::Max(i2-10,0),1);
3415         AliTPCseed * track2 = MakeSeed(track,0.2,0.5,0.9);
3416         FollowProlongation(*track2, i2,1);
3417         track2->SetBConstrain(kFALSE);
3418         track2->SetSeedType(4);
3419         arr->AddLast(track2);
3420         seed->Reset();
3421         seed->~AliTPCseed();
3422       }
3423
3424
3425       //arr->AddLast(track);
3426       //seed = new AliTPCseed;
3427       nout3++;
3428     }
3429   }
3430
3431   if (fDebug>3){
3432     Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
3433   }
3434   delete seed;
3435 }
3436
3437
3438 //_____________________________________________________________________________
3439 void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2, Float_t */*cuts[4]*/,
3440                                  Float_t deltay, Bool_t /*bconstrain*/) {
3441   //-----------------------------------------------------------------
3442   // This function creates track seeds - without vertex constraint
3443   //-----------------------------------------------------------------
3444   // cuts[0]   - fP4 cut        - not applied
3445   // cuts[1]   - tan(phi)  cut
3446   // cuts[2]   - zvertex cut    - not applied
3447   // cuts[3]   - fP3 cut
3448   Int_t nin0=0;
3449   Int_t nin1=0;
3450   Int_t nin2=0;
3451   Int_t nin3=0;
3452   //  Int_t nin4=0;
3453   //Int_t nin5=0;
3454
3455
3456
3457   Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
3458   //  Double_t cs=cos(alpha), sn=sin(alpha);
3459   Int_t row0 = (i1+i2)/2;
3460   Int_t drow = (i1-i2)/2;
3461   const AliTPCtrackerRow& kr0=fSectors[sec][row0];
3462   AliTPCtrackerRow * kr=0;
3463
3464   AliTPCpolyTrack polytrack;
3465   Int_t nclusters=fSectors[sec][row0];
3466   AliTPCseed * seed = new AliTPCseed;
3467
3468   Int_t sumused=0;
3469   Int_t cused=0;
3470   Int_t cnused=0;
3471   for (Int_t is=0; is < nclusters; is++) {  //LOOP over clusters
3472     Int_t nfound =0;
3473     Int_t nfoundable =0;
3474     for (Int_t iter =1; iter<2; iter++){   //iterations
3475       const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
3476       const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];
3477       const AliTPCclusterMI * cl= kr0[is];
3478
3479       if (cl->IsUsed(10)) {
3480         cused++;
3481       }
3482       else{
3483         cnused++;
3484       }
3485       Double_t x = kr0.GetX();
3486       // Initialization of the polytrack
3487       nfound =0;
3488       nfoundable =0;
3489       polytrack.Reset();
3490       //
3491       Double_t y0= cl->GetY();
3492       Double_t z0= cl->GetZ();
3493       Float_t erry = 0;
3494       Float_t errz = 0;
3495
3496       Double_t ymax = fSectors->GetMaxY(row0)-kr0.GetDeadZone()-1.5;
3497       if (deltay>0 && TMath::Abs(ymax-TMath::Abs(y0))> deltay ) continue;  // seed only at the edge
3498
3499       erry = (0.5)*cl->GetSigmaY2()/TMath::Sqrt(cl->GetQ())*6;
3500       errz = (0.5)*cl->GetSigmaZ2()/TMath::Sqrt(cl->GetQ())*6;
3501       polytrack.AddPoint(x,y0,z0,erry, errz);
3502
3503       sumused=0;
3504       if (cl->IsUsed(10)) sumused++;
3505
3506
3507       Float_t roady = (5*TMath::Sqrt(cl->GetSigmaY2()+0.2)+1.)*iter;
3508       Float_t roadz = (5*TMath::Sqrt(cl->GetSigmaZ2()+0.2)+1.)*iter;
3509       //
3510       x = krm.GetX();
3511       AliTPCclusterMI * cl1 = krm.FindNearest(y0,z0,roady,roadz);
3512       if (cl1 && TMath::Abs(ymax-TMath::Abs(y0))) {
3513         erry = (0.5)*cl1->GetSigmaY2()/TMath::Sqrt(cl1->GetQ())*3;
3514         errz = (0.5)*cl1->GetSigmaZ2()/TMath::Sqrt(cl1->GetQ())*3;
3515         if (cl1->IsUsed(10))  sumused++;
3516         polytrack.AddPoint(x,cl1->GetY(),cl1->GetZ(),erry,errz);
3517       }
3518       //
3519       x = krp.GetX();
3520       AliTPCclusterMI * cl2 = krp.FindNearest(y0,z0,roady,roadz);
3521       if (cl2) {
3522         erry = (0.5)*cl2->GetSigmaY2()/TMath::Sqrt(cl2->GetQ())*3;
3523         errz = (0.5)*cl2->GetSigmaZ2()/TMath::Sqrt(cl2->GetQ())*3;
3524         if (cl2->IsUsed(10)) sumused++;
3525         polytrack.AddPoint(x,cl2->GetY(),cl2->GetZ(),erry,errz);
3526       }
3527       //
3528       if (sumused>0) continue;
3529       nin0++;
3530       polytrack.UpdateParameters();
3531       // follow polytrack
3532       roadz = 1.2;
3533       roady = 1.2;
3534       //
3535       Double_t yn,zn;
3536       nfoundable = polytrack.GetN();
3537       nfound     = nfoundable;
3538       //
3539       for (Int_t ddrow = iter+1; ddrow<drow;ddrow++){
3540         Float_t maxdist = 0.8*(1.+3./(ddrow));
3541         for (Int_t delta = -1;delta<=1;delta+=2){
3542           Int_t row = row0+ddrow*delta;
3543           kr = &(fSectors[sec][row]);
3544           Double_t xn = kr->GetX();
3545           Double_t ymax = fSectors->GetMaxY(row)-kr->GetDeadZone()-1.5;
3546           polytrack.GetFitPoint(xn,yn,zn);
3547           if (TMath::Abs(yn)>ymax) continue;
3548           nfoundable++;
3549           AliTPCclusterMI * cln = kr->FindNearest(yn,zn,roady,roadz);
3550           if (cln) {
3551             Float_t dist =  TMath::Sqrt(  (yn-cln->GetY())*(yn-cln->GetY())+(zn-cln->GetZ())*(zn-cln->GetZ()));
3552             if (dist<maxdist){
3553               /*
3554               erry = (dist+0.3)*cln->GetSigmaY2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3555               errz = (dist+0.3)*cln->GetSigmaZ2()/TMath::Sqrt(cln->GetQ())*(1.+1./(ddrow));
3556               if (cln->IsUsed(10)) {
3557                 //      printf("used\n");
3558                 sumused++;
3559                 erry*=2;
3560                 errz*=2;
3561               }
3562               */
3563               erry=0.1;
3564               errz=0.1;
3565               polytrack.AddPoint(xn,cln->GetY(),cln->GetZ(),erry, errz);
3566               nfound++;
3567             }
3568           }
3569         }
3570         if ( (sumused>3) || (sumused>0.5*nfound) || (nfound<0.6*nfoundable))  break;
3571         polytrack.UpdateParameters();
3572       }
3573     }
3574     if ( (sumused>3) || (sumused>0.5*nfound))  {
3575       //printf("sumused   %d\n",sumused);
3576       continue;
3577     }
3578     nin1++;
3579     Double_t dy,dz;
3580     polytrack.GetFitDerivation(kr0.GetX(),dy,dz);
3581     AliTPCpolyTrack track2;
3582
3583     polytrack.Refit(track2,0.5+TMath::Abs(dy)*0.3,0.4+TMath::Abs(dz)*0.3);
3584     if (track2.GetN()<0.5*nfoundable) continue;