ITS/AliITStrackerSA.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////
  19 //  Stand alone ITS tracker class                        //
  20 //  Origin:  Elisabetta Crescio - crescio@to.infn.it     //
  21 //  Updated: Francesco Prino    - prino@to.infn.it       //
  22 ///////////////////////////////////////////////////////////
  23
  24 #include <stdlib.h>
  25
  26 #include <TArrayI.h>
  27 #include <TBranch.h>
  28 #include <TObjArray.h>
  29 #include <TTree.h>
  30
  31 #include "AliESDEvent.h"
  32 #include "AliESDVertex.h"
  33 #include "AliESDtrack.h"
  34 #include "AliITSVertexer.h"
  35 #include "AliITSclusterTable.h"
  36 #include "AliITSRecPoint.h"
  37 #include "AliITSgeomTGeo.h"
  38 #include "AliITStrackSA.h"
  39 #include "AliITStrackerSA.h"
  40 #include "AliITSReconstructor.h"
  41 #include "AliLog.h"
  42 #include "AliRun.h"
  43
  44 ClassImp(AliITStrackerSA)
  45
  46 //____________________________________________________________________________
  47 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(),
  48 fPhiEstimate(0),
  49 fITSStandAlone(0),
  50 fLambdac(0),
  51 fPhic(0),
  52 fCoef1(0),
  53 fCoef2(0),
  54 fCoef3(0),
  55 fNloop(0),
  56 fPhiWin(0),
  57 fLambdaWin(0),
  58 fListOfTracks(0),
  59 fListOfSATracks(0),
  60 fITSclusters(0),
  61 fInwardFlag(0),
  62 fOuterStartLayer(0),
  63 fInnerStartLayer(5),
  64 fMinNPoints(0),
  65 fMinQ(0.),
  66 fCluCoord(0){
  67   // Default constructor
  68   Init();
  69
  70 }
  71 //____________________________________________________________________________
  72 AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
  73 fPhiEstimate(0),
  74 fITSStandAlone(0),
  75 fLambdac(0),
  76 fPhic(0),
  77 fCoef1(0),
  78 fCoef2(0),
  79 fCoef3(0),
  80 fNloop(0),
  81 fPhiWin(0),
  82 fLambdaWin(0),
  83 fListOfTracks(0),
  84 fListOfSATracks(0),
  85 fITSclusters(0),
  86 fInwardFlag(0),
  87 fOuterStartLayer(0),
  88 fInnerStartLayer(5),
  89 fMinNPoints(0),
  90 fMinQ(0.),
  91 fCluCoord(0)
  92 {
  93   // Standard constructor (Vertex is known and passed to this obj.)
  94   if (geom) {
  95     AliWarning("\"geom\" is actually a dummy argument !");
  96   }
  97
  98   Init();
  99
 100 }
 101
 102 //____________________________________________________________________________
 103 AliITStrackerSA::~AliITStrackerSA(){
 104   // destructor
 105
 106   if(fPhiWin)delete []fPhiWin;
 107   if(fLambdaWin)delete []fLambdaWin;
 108   fListOfTracks->Delete();
 109   delete fListOfTracks;
 110   fListOfSATracks->Delete();
 111   delete fListOfSATracks;
 112   if(fCluCoord){
 113     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 114       if(fCluCoord[i]){
 115         fCluCoord[i]->Delete();
 116         delete fCluCoord[i];
 117       }
 118     }
 119     delete [] fCluCoord;
 120   }
 121 }
 122
 123 //____________________________________________________________________________
 124 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
 125 // This method is used to find and fit the tracks. By default the corresponding
 126 // method in the parent class is invoked. In this way a combined tracking
 127 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
 128 // is done in the ITS only. In the standard reconstruction chain this option
 129 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
 130   Int_t rc=0;
 131
 132   if(!fITSStandAlone){
 133     rc=AliITStrackerMI::Clusters2Tracks(event);
 134   }
 135   else {
 136     AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
 137   }
 138   if(!rc){
 139     rc=FindTracks(event,kFALSE);
 140     Int_t nSPDcontr=0;
 141     const AliESDVertex *spdv = event->GetPrimaryVertexSPD();
 142     if(spdv) nSPDcontr = spdv->GetNContributors();
 143     if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE &&
 144        nSPDcontr<=AliITSReconstructor::GetRecoParam()->GetMaxSPDcontrForSAToUseAllClusters()) {
 145       rc=FindTracks(event,kTRUE);
 146     }
 147   }
 148   return rc;
 149 }
 150
 151 //____________________________________________________________________________
 152 void AliITStrackerSA::Init(){
 153   //  Reset all data members
 154     fPhiEstimate=0;
 155     for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 156     fLambdac=0;
 157     fPhic=0;
 158     fCoef1=0;
 159     fCoef2=0;
 160     fCoef3=0;
 161     fPointc[0]=0;
 162     fPointc[1]=0;
 163     Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
 164     if(nLoops==32){
 165       SetFixedWindowSizes();
 166     }else{
 167       Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
 168       Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
 169       Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
 170       Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
 171       SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
 172     }
 173     fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
 174     fITSclusters = 0;
 175     SetOuterStartLayer(1);
 176     SetSAFlag(kFALSE);
 177     fListOfTracks=new TClonesArray("AliITStrackMI",100);
 178     fListOfSATracks=new TClonesArray("AliITStrackSA",100);
 179     fCluCoord = 0;
 180     fMinNPoints = 3;
 181  }
 182 //_______________________________________________________________________
 183 void AliITStrackerSA::ResetForFinding(){
 184   //  Reset data members used in all loops during track finding
 185     fPhiEstimate=0;
 186     for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 187     fLambdac=0;
 188     fPhic=0;
 189     fCoef1=0;
 190     fCoef2=0;
 191     fCoef3=0;
 192     fPointc[0]=0;
 193     fPointc[1]=0;
 194     fListOfTracks->Clear();
 195     fListOfSATracks->Clear();
 196 }
 197
 198
 199
 200 //______________________________________________________________________
 201 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event, Bool_t useAllClusters){
 202
 203 // Track finder using the ESD object
 204
 205   AliDebug(2,Form(" field is %f",event->GetMagneticField()));
 206   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
 207
 208   if(!fITSclusters){
 209     Fatal("FindTracks","ITS cluster tree is not accessed!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
 210     return -1;
 211   }
 212   //Reads event and mark clusters of traks already found, with flag kITSin
 213   Int_t nentr=event->GetNumberOfTracks();
 214   if(!useAllClusters) {
 215     while (nentr--) {
 216       AliESDtrack *track=event->GetTrack(nentr);
 217       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
 218         Int_t idx[12];
 219         Int_t ncl = track->GetITSclusters(idx);
 220         for(Int_t k=0;k<ncl;k++){
 221           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
 222           cll->SetBit(kSAflag);
 223         }
 224       }
 225     }
 226   }else{
 227     while (nentr--) {
 228       AliESDtrack *track=event->GetTrack(nentr);
 229       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
 230         Int_t idx[12];
 231         Int_t ncl = track->GetITSclusters(idx);
 232         for(Int_t k=0;k<ncl;k++){
 233           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
 234           cll->ResetBit(kSAflag);
 235         }
 236       }
 237     }
 238   }
 239   //Get primary vertex
 240   Double_t primaryVertex[3];
 241   event->GetVertex()->GetXYZ(primaryVertex);
 242   //Creates TClonesArray with clusters for each layer. The clusters already used
 243   //by AliITStrackerMI are not considered
 244   Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
 245   Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
 246   if (fCluCoord == 0) {
 247     fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
 248     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 249       fCluCoord[i]=0;
 250     }
 251   }
 252   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 253     AliITSlayer &layer=fgLayers[i];
 254     if (!ForceSkippingOfLayer(i)) {
 255       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 256         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 257         if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
 258         if(cls->GetQ()==0) continue; //fake clusters dead zones
 259         if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
 260         nclusters[i]++;
 261       }
 262     }
 263     dmar[i]=0;
 264     if(!fCluCoord[i]){
 265       fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
 266     }else{
 267       fCluCoord[i]->Delete();
 268       fCluCoord[i]->Expand(nclusters[i]);
 269     }
 270   }
 271
 272   for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
 273     TClonesArray &clucoo = *fCluCoord[ilay];
 274     AliITSlayer &layer=fgLayers[ilay];
 275     if (!ForceSkippingOfLayer(ilay)) {
 276       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 277         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 278         if(cls->TestBit(kSAflag)==kTRUE) continue;
 279         if(cls->GetQ()==0) continue;
 280         if(ilay>1 && cls->GetQ()<=fMinQ) continue;
 281         Double_t phi=0;Double_t lambda=0;
 282         Double_t x=0;Double_t y=0;Double_t z=0;
 283         Double_t sx=0;Double_t sy=0;Double_t sz=0;
 284         GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
 285         GetCoorErrors(cls,sx,sy,sz);
 286         new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
 287         dmar[ilay]++;
 288       }
 289     }
 290     fCluCoord[ilay]->Sort();
 291   }
 292
 293   // track counter
 294   Int_t ntrack=0;
 295
 296   static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
 297   Int_t startLayForSeed=0;
 298   Int_t lastLayForSeed=fOuterStartLayer;
 299   Int_t nSeedSteps=lastLayForSeed-startLayForSeed;
 300   Int_t seedStep=1;
 301   if(fInwardFlag){
 302     startLayForSeed=AliITSgeomTGeo::GetNLayers()-1;
 303     lastLayForSeed=fInnerStartLayer;
 304     nSeedSteps=startLayForSeed-lastLayForSeed;
 305     seedStep=-1;
 306   }
 307
 308   // loop on minimum number of points
 309   for(Int_t iMinNPoints=AliITSgeomTGeo::GetNLayers(); iMinNPoints>=fMinNPoints; iMinNPoints--) {
 310
 311     // loop on starting layer for track finding
 312     for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) {
 313       Int_t theLay=startLayForSeed+iSeedLay*seedStep;
 314       if(ForceSkippingOfLayer(theLay)) continue;
 315       Int_t minNPoints=iMinNPoints-theLay;
 316       if(fInwardFlag) minNPoints=iMinNPoints-(AliITSgeomTGeo::GetNLayers()-1-theLay);
 317       for(Int_t i=theLay+1;i<AliITSgeomTGeo::GetNLayers();i++)
 318         if(ForceSkippingOfLayer(i))
 319           minNPoints--;
 320       if(minNPoints<fMinNPoints) continue;
 321
 322       // loop on phi and lambda window size
 323       for(Int_t nloop=0;nloop<fNloop;nloop++){
 324         Int_t nclTheLay=fCluCoord[theLay]->GetEntries();
 325         while(nclTheLay--){
 326           ResetForFinding();
 327           Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex);
 328           if(!useRP) continue;
 329           AliITStrackSA trs;
 330
 331           Int_t pflag=0;
 332           Int_t kk;
 333           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
 334
 335           kk=0;
 336           nClusLay[kk] = SearchClusters(theLay,fPhiWin[nloop],fLambdaWin[nloop],
 337                                         &trs,primaryVertex[2],pflag);
 338           Int_t nextLay=theLay+seedStep;
 339           Bool_t goon=kTRUE;
 340           if(nextLay<0 || nextLay == 6) goon = kFALSE;
 341           while(goon){
 342             kk++;
 343             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 344                                             &trs,primaryVertex[2],pflag);
 345             if(nClusLay[kk]!=0){
 346               pflag=1;
 347               if(kk==1) {
 348                 fPoint3[0]=fPointc[0];
 349                 fPoint3[1]=fPointc[1];
 350               } else {
 351                 UpdatePoints();
 352               }
 353             }
 354             nextLay+=seedStep;
 355             if(nextLay<0 || nextLay==6) goon=kFALSE;
 356           }
 357
 358
 359           Int_t layOK=0;
 360           if(!fInwardFlag){
 361             for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-theLay;nnp++){
 362               if(nClusLay[nnp]!=0) layOK+=1;
 363             }
 364           }else{
 365             for(Int_t nnp=theLay; nnp>=0; nnp--){
 366               if(nClusLay[nnp]!=0) layOK+=1;
 367             }
 368           }
 369           if(layOK>=minNPoints){
 370             AliDebug(2,Form("---NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]));
 371             AliITStrackV2* tr2 = 0;
 372             tr2 = FitTrack(&trs,primaryVertex);
 373             if(!tr2){
 374               continue;
 375             }
 376             AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
 377
 378             StoreTrack(tr2,event,useAllClusters);
 379             ntrack++;
 380
 381           }
 382
 383         }//end loop on clusters of theLay
 384       } //end loop on window sizes
 385     } //end loop on theLay
 386   }//end loop on min points
 387
 388   // search for 1-point tracks in SPD, only for cosmics
 389   // (A.Dainese 21.03.08)
 390   if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
 391      TMath::Abs(event->GetMagneticField())<0.01) {
 392     Int_t outerLayer=1; // only SPD
 393     for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
 394       //   counter for clusters on each layer
 395
 396       for(Int_t nloop=0;nloop<fNloop;nloop++){
 397         Int_t nclInnLay=fCluCoord[innLay]->GetEntries();
 398         while(nclInnLay--){ //loop starting from layer innLay
 399           ResetForFinding();
 400           Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
 401           if(!useRP) continue;
 402           AliITStrackSA trs;
 403
 404           Int_t pflag=0;
 405           Int_t kk;
 406           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
 407
 408           kk=0;
 409           nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
 410                                   &trs,primaryVertex[2],pflag);
 411           for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
 412             kk++;
 413             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 414                                     &trs,primaryVertex[2],pflag);
 415             if(nClusLay[kk]!=0){
 416               pflag=1;
 417               if(kk==1) {
 418                 fPoint3[0]=fPointc[0];
 419                 fPoint3[1]=fPointc[1];
 420               } else {
 421                 UpdatePoints();
 422               }
 423             }
 424           }
 425
 426           Int_t layOK=0;
 427           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
 428             if(nClusLay[nnp]!=0) layOK+=1;
 429           }
 430           if(layOK==1) {
 431             AliDebug(2,Form("----NPOINTS: %d; MAP: %d %d %d %d %d %d\n",layOK,nClusLay[0],nClusLay[1],nClusLay[2],nClusLay[3],nClusLay[4],nClusLay[5]));
 432             AliITStrackV2* tr2 = 0;
 433             Bool_t onePoint = kTRUE;
 434             tr2 = FitTrack(&trs,primaryVertex,onePoint);
 435             if(!tr2){
 436               continue;
 437             }
 438             AliDebug(2,Form("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
 439
 440             StoreTrack(tr2,event,useAllClusters);
 441             ntrack++;
 442
 443           }
 444
 445         }//end loop on clusters of innLay
 446       } //end loop on window sizes
 447
 448     } //end loop on innLay
 449   } // end search 1-point tracks
 450
 451   if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
 452   ResetForFinding();
 453   return 0;
 454
 455 }
 456
 457 //________________________________________________________________________
 458
 459 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
 460   //fit of the found track (most general case, also <6 points, layers missing)
 461   // A.Dainese 16.11.07
 462
 463
 464   const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
 465
 466   static Int_t firstmod[AliITSgeomTGeo::kNLayers];
 467   static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
 468   static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
 469   static Int_t end[AliITSgeomTGeo::kNLayers];
 470   static Int_t indices[AliITSgeomTGeo::kNLayers];
 471
 472   static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
 473
 474   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 475     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
 476     end[i]=0;
 477     for(Int_t j=0;j<kMaxClu; j++){
 478       clind[i][j]=0;
 479       clmark[i][j]=0;
 480       listlayer[i][j]=0;
 481    }
 482   }
 483
 484
 485   Int_t nclusters = tr->GetNumberOfClustersSA();
 486   for(Int_t ncl=0;ncl<nclusters;ncl++){
 487     Int_t index = tr->GetClusterIndexSA(ncl);
 488     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
 489     Int_t lay = (index & 0xf0000000) >> 28;
 490     Int_t nInLay=end[lay];
 491     listlayer[lay][nInLay]=cl;
 492     clind[lay][nInLay]=index;
 493     end[lay]++;
 494   }
 495
 496   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
 497     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
 498       Int_t mark = tr->GetClusterMark(nlay,ncl);
 499       clmark[nlay][ncl]=mark;
 500     }
 501   }
 502
 503
 504   Int_t firstLay=-1,secondLay=-1;
 505   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 506     if(end[i]==0) {
 507       end[i]=1;
 508     }else{
 509       if(firstLay==-1) {
 510         firstLay=i;
 511       } else if(secondLay==-1) {
 512         secondLay=i;
 513       }
 514     }
 515   }
 516
 517   if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
 518   TClonesArray &arrMI= *fListOfTracks;
 519   TClonesArray &arrSA= *fListOfSATracks;
 520   Int_t nFoundTracks=0;
 521
 522
 523   for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
 524     indices[0]=l0;
 525     for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
 526       indices[1]=l1;
 527       for(Int_t l2=0;l2<end[2];l2++){  //loop on layer 3
 528         indices[2]=l2;
 529         for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
 530           indices[3]=l3;
 531           for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
 532             indices[4]=l4;
 533             for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
 534               indices[5]=l5;
 535
 536               // estimate curvature from 2 innermost points (or innermost point + vertex)
 537
 538               Int_t iFirstLay=indices[firstLay];
 539               Int_t mrk1=clmark[firstLay][iFirstLay];
 540
 541               AliITSRecPoint* p1=(AliITSRecPoint*)listlayer[firstLay][iFirstLay];
 542               Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
 543               Int_t layer,ladder,detector;
 544               AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
 545               Double_t yclu1 = p1->GetY();
 546               Double_t zclu1 = p1->GetZ();
 547
 548               Double_t x1,y1,z1;
 549               Double_t x2,y2,z2;
 550               Double_t cv=0,tgl2=0,phi2=0;
 551               AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1);
 552               x1 = arr1->GetX();
 553               y1 = arr1->GetY();
 554               z1 = arr1->GetZ();
 555
 556               if(secondLay>0) {
 557                 Int_t iSecondLay=indices[secondLay];
 558                 Int_t mrk2=clmark[secondLay][iSecondLay];
 559                 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2);
 560                 x2 = arr2->GetX();
 561                 y2 = arr2->GetY();
 562                 z2 = arr2->GetZ();
 563                 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
 564                 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
 565                 phi2 = TMath::ATan2((y2-y1),(x2-x1));
 566               } else { // special case of 1-point tracks, only for cosmics (B=0)
 567                 x2 = primaryVertex[0];
 568                 y2 = primaryVertex[1];
 569                 z2 = primaryVertex[2];
 570                 cv = 0;
 571                 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 572                 phi2 = TMath::ATan2((y1-y2),(x1-x2));
 573               }
 574
 575               // create track and attach it the RecPoints
 576               AliITStrackSA trac(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
 577               for(Int_t iLay=5; iLay>=0; iLay--){
 578                 Int_t iInLay=indices[iLay];
 579                 AliITSRecPoint* cl=(AliITSRecPoint*)listlayer[iLay][iInLay];
 580                 if(cl!=0){
 581                   trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0);
 582                   trac.AddClusterMark(iLay,clmark[iLay][iInLay]);
 583                 }
 584               }
 585
 586               //fit with Kalman filter using AliITStrackerMI::RefitAt()
 587               AliITStrackSA ot(trac);
 588
 589               ot.ResetCovariance(10.);
 590               ot.ResetClusters();
 591
 592               // Propagate inside the innermost layer with a cluster
 593               if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
 594
 595                 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&ot,&trac)){ //fit from layer 1 to layer 6
 596                   AliITStrackMI otrack2(ot);
 597                   otrack2.ResetCovariance(10.);
 598                   otrack2.ResetClusters();
 599                   //fit from layer 6 to layer 1
 600                   if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),&otrack2,&ot)) {
 601                     new(arrMI[nFoundTracks]) AliITStrackMI(otrack2);
 602                     new(arrSA[nFoundTracks]) AliITStrackSA(trac);
 603                     ++nFoundTracks;
 604                   }
 605
 606                 }
 607               }
 608             }//end loop layer 6
 609           }//end loop layer 5
 610         }//end loop layer 4
 611       }//end loop layer 3
 612     }//end loop layer 2
 613   }//end loop layer 1
 614
 615
 616
 617
 618   if(fListOfTracks->GetEntries()==0) return 0;
 619
 620   Int_t lowchi2 = FindTrackLowChiSquare();
 621   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
 622   AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
 623
 624   if(otrack==0) return 0;
 625
 626   CookLabel(otrack,0.); //MI change - to see fake ratio
 627   Int_t label=FindLabel(otrack);
 628   otrack->SetLabel(label);
 629   Double_t low=0.;
 630   Double_t up=0.51;
 631   otrack->CookdEdx(low,up);
 632
 633   //remove clusters of found track
 634   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
 635     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
 636       Int_t index = trsa->GetClusterMark(nlay,cln);
 637       RemoveClusterCoord(nlay,index);
 638     }
 639   }
 640
 641   return otrack;
 642
 643 }
 644
 645 //_______________________________________________________
 646 void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) const
 647 {
 648   //
 649   // Add new track to the ESD
 650   //
 651   AliESDtrack outtrack;
 652   outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
 653   if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
 654   for(Int_t i=0;i<12;i++) {
 655     outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
 656   }
 657   Double_t sdedx[4]={0.,0.,0.,0.};
 658   for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
 659   outtrack.SetITSdEdxSamples(sdedx);
 660
 661
 662   if(AliITSReconstructor::GetRecoParam()->GetSAUsedEdxInfo()){
 663     Double_t mom=t->P();
 664     Double_t ppid[AliPID::kSPECIES];
 665     for(Int_t isp=0;isp<AliPID::kSPECIES;isp++) ppid[isp]=0.;
 666     ppid[AliPID::kPion]=1.;
 667     if(mom<0.7){
 668       Double_t truncmean=t->GetdEdx();
 669       Int_t ide=fITSPid->GetParticleIdFromdEdxVsP(mom,truncmean,kTRUE);
 670       if(ide==AliPID::kProton){
 671         ppid[AliPID::kProton]=1.;
 672         ppid[AliPID::kPion]=0.;
 673       }
 674       else if(ide==AliPID::kKaon){
 675         ppid[AliPID::kKaon]=1.;
 676         ppid[AliPID::kPion]=0.;
 677       }
 678     }
 679     outtrack.SetITSpid(ppid);
 680     outtrack.SetESDpid(ppid);
 681   }
 682   event->AddTrack(&outtrack);
 683
 684   return;
 685 }
 686
 687
 688 //_______________________________________________________
 689 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
 690   //function used to to find the clusters associated to the track
 691
 692   if(ForceSkippingOfLayer(layer)) return 0;
 693
 694
 695   Int_t nc=0;
 696   AliITSlayer &lay = fgLayers[layer];
 697   Double_t r=lay.GetR();
 698   if(pflag==1){
 699     Double_t cx1,cx2,cy1,cy2;
 700     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
 701     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
 702        return 0;
 703     Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
 704     Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
 705     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
 706   }
 707
 708
 709   Double_t phiExpect=fPhiEstimate;
 710   Double_t lamExpect=fLambdac;
 711
 712   Int_t ncl = fCluCoord[layer]->GetEntriesFast();
 713   Int_t startcl=FindIndex(layer,lamExpect-lambdawindow*1.02);
 714   Int_t endcl=FindIndex(layer,lamExpect+lambdawindow*1.02)+1;
 715   if(endcl>=ncl) endcl=ncl-1;
 716
 717   for (Int_t index=startcl; index<=endcl; index++) {
 718     //for (Int_t index=0; index<ncl; index++) {
 719     AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
 720
 721
 722     Double_t phi = arr->GetPhi();
 723     Double_t deltaPhi = phi-phiExpect;
 724     if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi();
 725     else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi();
 726     if (TMath::Abs(deltaPhi)>phiwindow) continue;
 727
 728     Double_t lambda = arr->GetLambda();
 729     if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue;
 730
 731     if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
 732     if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
 733     Int_t orind = arr->GetOrInd();
 734     trs->AddClusterSA(layer,orind);
 735     trs->AddClusterMark(layer,index);
 736     nc++;
 737     fLambdac=lambda;
 738     fPhiEstimate=phi;
 739
 740     fPointc[0]=arr->GetX();
 741     fPointc[1]=arr->GetY();
 742
 743   }
 744   return nc;
 745 }
 746
 747 //________________________________________________________________
 748 Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
 749   // Sets the first point (seed) for tracking
 750
 751   AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
 752   fPhic = arr->GetPhi();
 753   fLambdac = arr->GetLambda();
 754   fPhiEstimate = fPhic;
 755   fPoint1[0]=primaryVertex[0];
 756   fPoint1[1]=primaryVertex[1];
 757   fPoint2[0]=arr->GetX();
 758   fPoint2[1]=arr->GetY();
 759   return kTRUE;
 760 }
 761
 762 //________________________________________________________________
 763 void AliITStrackerSA::UpdatePoints(){
 764   //update of points for the estimation of the curvature
 765
 766   fPoint2[0]=fPoint3[0];
 767   fPoint2[1]=fPoint3[1];
 768   fPoint3[0]=fPointc[0];
 769   fPoint3[1]=fPointc[1];
 770
 771
 772 }
 773
 774 //___________________________________________________________________
 775 Int_t AliITStrackerSA::FindEquation(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3,Double_t& a, Double_t& b, Double_t& c){
 776
 777    //given (x,y) of three recpoints (in global coordinates)
 778    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
 779    double dx31=x3-x1,dy31=y3-y1,dx21=x2-x1,dy21=y2-y1;
 780    Double_t den = dx31*dy21-dx21*dy31;
 781    if(den==0) return 0;
 782    else den = 1./den;
 783    //
 784    double r31 = -dx31*(x1+x3) - dy31*(y1+y3);
 785    double r21 = -dx21*(x1+x2) - dy21*(y1+y2);
 786    //
 787    double da = r31*dy21 - r21*dy31;
 788    double db = r21*dx31 - r31*dx21;
 789    a = da*den;
 790    b = db*den;
 791    c = -x1*x1-y1*y1-a*x1-b*y1;
 792    return 1;
 793  }
 794 //__________________________________________________________________________
 795  Int_t AliITStrackerSA::FindIntersection(Double_t a1, Double_t b1, Double_t c1, Double_t c2,Double_t& x1,Double_t& y1, Double_t& x2, Double_t& y2){
 796
 797  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
 798  //c2 is -rlayer*rlayer
 799
 800   if(a1==0) return 0;
 801  Double_t m = c2-c1;
 802  Double_t aA = (b1*b1)/(a1*a1)+1;
 803  Double_t bB = (-2*m*b1/(a1*a1));
 804  Double_t cC = c2+(m*m)/(a1*a1);
 805  Double_t dD = bB*bB-4*aA*cC;
 806  if(dD<0) return 0;
 807
 808  y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
 809  y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
 810  x1 = (c2-c1-b1*y1)/a1;
 811  x2 = (c2-c1-b1*y2)/a1;
 812
 813  return 1;
 814 }
 815 //____________________________________________________________________
 816 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
 817 x2,Double_t y2,Double_t x3,Double_t y3){
 818
 819   //calculates the curvature of track
 820   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 821   if(den==0) return 0;
 822   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 823   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 824   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
 825   Double_t xc=-a/2.;
 826
 827   if((a*a+b*b-4*c)<0) return 0;
 828   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
 829   if(rad==0) return 0;
 830
 831   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
 832   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
 833   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
 834   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
 835
 836   return 1/rad;
 837
 838 }
 839
 840
 841 //____________________________________________________________________
 842 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
 843
 844   //Returns the point closest to pp
 845
 846   Double_t diff1 = p1-pp;
 847   Double_t diff2 = p2-pp;
 848
 849   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
 850   else fPhiEstimate=p2;
 851   return fPhiEstimate;
 852
 853 }
 854
 855
 856 //_________________________________________________________________
 857 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
 858   // returns track with lowest chi square
 859   Int_t dim=fListOfTracks->GetEntries();
 860   if(dim<=1) return 0;
 861   AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
 862   Double_t minChi2=trk->GetChi2();
 863   Int_t index=0;
 864   for(Int_t i=1;i<dim;i++){
 865     trk = (AliITStrackV2*)fListOfTracks->At(i);
 866     Double_t chi2=trk->GetChi2();
 867     if(chi2<minChi2){
 868       minChi2=chi2;
 869       index=i;
 870     }
 871   }
 872   return index;
 873 }
 874
 875 //__________________________________________________________
 876 Int_t AliITStrackerSA::FindLabel(AliITStrackV2* track){
 877   // compute the track label starting from cluster labels
 878
 879   Int_t labl[AliITSgeomTGeo::kNLayers][3];
 880   Int_t cnts[AliITSgeomTGeo::kNLayers][3];
 881   for(Int_t j=0;j<AliITSgeomTGeo::GetNLayers();j++){
 882     for(Int_t k=0;k<3;k++){
 883       labl[j][k]=-2;
 884       cnts[j][k]=1;
 885     }
 886   }
 887   Int_t iNotLabel=0;
 888   for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
 889     Int_t indexc = track->GetClusterIndex(i);
 890     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc);
 891     Int_t iLayer=cl->GetLayer();
 892     for(Int_t k=0;k<3;k++){
 893       labl[iLayer][k]=cl->GetLabel(k);
 894       if(labl[iLayer][k]<0) iNotLabel++;
 895     }
 896   }
 897   if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
 898
 899   for(Int_t j1=0;j1<AliITSgeomTGeo::kNLayers; j1++) {
 900     for(Int_t j2=0; j2<j1;  j2++){
 901       for(Int_t k1=0; k1<3; k1++){
 902         for(Int_t k2=0; k2<3; k2++){
 903           if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
 904             cnts[j2][k2]++;
 905             cnts[j1][k1]=0;
 906           }
 907         }
 908       }
 909     }
 910   }
 911
 912
 913   Int_t cntMax=0;
 914   Int_t label=-1;
 915   for(Int_t j=0;j<AliITSgeomTGeo::kNLayers;j++){
 916     for(Int_t k=0;k<3;k++){
 917       if(cnts[j][k]>cntMax && labl[j][k]>=0){
 918         cntMax=cnts[j][k];
 919         label=labl[j][k];
 920       }
 921     }
 922   }
 923
 924   Int_t lflag=0;
 925   for(Int_t i=0;i<AliITSgeomTGeo::kNLayers;i++)
 926     if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
 927
 928   if(lflag<track->GetNumberOfClusters()) label = -label;
 929   return label;
 930 }
 931 //_____________________________________________________________________________
 932 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Double_t phimin, Double_t phimax, Double_t lambdamin, Double_t lambdamax){
 933   // Set sizes of the phi and lambda windows used for track finding
 934   fNloop = n;
 935   if(fPhiWin) delete [] fPhiWin;
 936   if(fLambdaWin) delete [] fLambdaWin;
 937   fPhiWin = new Double_t[fNloop];
 938   fLambdaWin = new Double_t[fNloop];
 939   Double_t stepPhi=(phimax-phimin)/(Double_t)(fNloop-1);
 940   Double_t stepLambda=(lambdamax-lambdamin)/(Double_t)(fNloop-1);
 941   for(Int_t k=0;k<fNloop;k++){
 942     Double_t phi=phimin+k*stepPhi;
 943     Double_t lam=lambdamin+k*stepLambda;
 944     fPhiWin[k]=phi;
 945     fLambdaWin[k]=lam;
 946   }
 947 }
 948 //_____________________________________________________________________________
 949 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
 950   // Set sizes of the phi and lambda windows used for track finding
 951   fNloop = n;
 952   if(phi){ // user defined values
 953     fPhiWin = new Double_t[fNloop];
 954     fLambdaWin = new Double_t[fNloop];
 955     for(Int_t k=0;k<fNloop;k++){
 956       fPhiWin[k]=phi[k];
 957       fLambdaWin[k]=lam[k];
 958     }
 959   }
 960   else {  // default values
 961
 962     Double_t phid[32]   = {0.002,0.003,0.004,0.0045,0.0047,
 963                            0.005,0.0053,0.0055,0.006,0.0063,
 964                            0.0065,0.007,0.0073,0.0075,0.0077,
 965                            0.008,0.0083,0.0085,0.0087,0.009,
 966                            0.0095,0.0097,0.01,0.0105,0.011,
 967                            0.0115,0.012,0.0125,0.013,0.0135,
 968                            0.0140,0.0145};
 969     Double_t lambdad[32] = {0.003,0.004,0.005,0.005,0.005,
 970                             0.005,0.005,0.006,0.006,0.006,
 971                             0.006,0.007,0.007,0.007,0.007,
 972                             0.007,0.007,0.007,0.007,0.007,
 973                             0.007,0.007,0.008,0.008,0.008,
 974                             0.008,0.008,0.008,0.008,0.008,
 975                             0.008,0.008};
 976
 977     if(fNloop!=32){
 978       fNloop = 32;
 979     }
 980
 981
 982     fPhiWin = new Double_t[fNloop];
 983     fLambdaWin = new Double_t[fNloop];
 984
 985     Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
 986
 987     for(Int_t k=0;k<fNloop;k++){
 988       fPhiWin[k]=phid[k]*factor;
 989       fLambdaWin[k]=lambdad[k]*factor;
 990     }
 991
 992   }
 993
 994 }
 995 //_______________________________________________________________________
 996 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Double_t &x, Double_t &y,Double_t &z, const Double_t* vertex){
 997   //Returns values of phi (azimuthal) and lambda angles for a given cluster
 998 /*
 999   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1000   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1001   Double_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1002
1003   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1004   Double_t phi1=TMath::Pi()/2+alpha;
1005   if (lay==1) phi1+=TMath::Pi();
1006
1007   Double_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1008   Double_t r=tx*cp+ty*sp;
1009
1010   xyz= r*cp - cl->GetY()*sp;
1011   y= r*sp + cl->GetY()*cp;
1012   z=cl->GetZ();
1013 */
1014   Float_t xyz[3];
1015   cl->GetGlobalXYZ(xyz);
1016   x=xyz[0];
1017   y=xyz[1];
1018   z=xyz[2];
1019
1020   phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1021   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1022 }
1023
1024 //________________________________________________________________________
1025 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Double_t &sx,Double_t &sy, Double_t &sz){
1026
1027   //returns sigmax, y, z of cluster in global coordinates
1028 /*
1029   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1030   Int_t lay,lad,det;
1031   AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1032
1033   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1034   Double_t phi=TMath::Pi()/2+alpha;
1035   if (lay==1) phi+=TMath::Pi();
1036
1037   Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1038 */
1039   Float_t covm[6];
1040   cl->GetGlobalCov(covm);
1041   sx=TMath::Sqrt(covm[0]);
1042   sy=TMath::Sqrt(covm[3]);
1043   sz=TMath::Sqrt(covm[5]);
1044 /*
1045   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1046   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1047   sz = TMath::Sqrt(cl->GetSigmaZ2());
1048 */
1049 }
1050
1051 //________________________________________________________________________
1052 Int_t AliITStrackerSA::FindIndex(Int_t lay, Double_t lamVal) const {
1053   // Find the cluster at limit of lambda window
1054
1055   Int_t base = 0;
1056   Int_t last = fCluCoord[lay]->GetEntriesFast()-1;
1057   if(last<0) return 0;
1058   Int_t position;
1059   Double_t lamfirst=((AliITSclusterTable*)fCluCoord[lay]->At(base))->GetLambda();
1060   if(lamfirst>lamVal) return base;
1061   Double_t lamlast=((AliITSclusterTable*)fCluCoord[lay]->At(last))->GetLambda();
1062   if(lamlast<=lamVal) return last;
1063   while (last >= base) {
1064     position = (base+last) / 2;
1065     Double_t a=((AliITSclusterTable*)fCluCoord[lay]->At(position))->GetLambda()-lamVal;
1066     Double_t b=((AliITSclusterTable*)fCluCoord[lay]->At(position+1))->GetLambda()-lamVal;
1067     if(a*b<=0) return position;
1068     if(a>0) last = position;
1069     else  base = position;
1070   }
1071   return 0;
1072 }
1073