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 tracker class                     //
  20 //  Origin:  Elisabetta Crescio                   //
  21 //  e-mail:  crescio@to.infn.it                   //
  22 //  tracks are saved as AliITStrackV2 objects     //
  23 ////////////////////////////////////////////////////
  24
  25 #include <stdlib.h>
  26
  27 #include <TArrayI.h>
  28 #include <TBranch.h>
  29 #include <TObjArray.h>
  30 #include <TTree.h>
  31
  32 #include "AliESDEvent.h"
  33 #include "AliESDVertex.h"
  34 #include "AliESDtrack.h"
  35 #include "AliITSVertexer.h"
  36 #include "AliITSclusterTable.h"
  37 #include "AliITSRecPoint.h"
  38 #include "AliITSgeomTGeo.h"
  39 #include "AliITStrackSA.h"
  40 #include "AliITStrackerSA.h"
  41 #include "AliITSReconstructor.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 fVert(0),
  59 fVertexer(0),
  60 fListOfTracks(0),
  61 fITSclusters(0),
  62 fSixPoints(0),
  63 fOuterStartLayer(0),
  64 fCluLayer(0),
  65 fCluCoord(0){
  66   // Default constructor
  67   Init();
  68
  69 }
  70 //____________________________________________________________________________
  71 AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
  72 fPhiEstimate(0),
  73 fITSStandAlone(0),
  74 fLambdac(0),
  75 fPhic(0),
  76 fCoef1(0),
  77 fCoef2(0),
  78 fCoef3(0),
  79 fNloop(0),
  80 fPhiWin(0),
  81 fLambdaWin(0),
  82 fVert(0),
  83 fVertexer(0),
  84 fListOfTracks(0),
  85 fITSclusters(0),
  86 fSixPoints(0),
  87 fOuterStartLayer(0),
  88 fCluLayer(0),
  89 fCluCoord(0)
  90 {
  91   // Standard constructor (Vertex is known and passed to this obj.)
  92   if (geom) {
  93     AliWarning("\"geom\" is actually a dummy argument !");
  94   }
  95
  96   Init();
  97   fVert = 0;
  98
  99 }
 100
 101 //____________________________________________________________________________
 102 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliESDVertex *vert):AliITStrackerMI(0),
 103 fPhiEstimate(0),
 104 fITSStandAlone(0),
 105 fLambdac(0),
 106 fPhic(0),
 107 fCoef1(0),
 108 fCoef2(0),
 109 fCoef3(0),
 110 fNloop(0),
 111 fPhiWin(0),
 112 fLambdaWin(0),
 113 fVert(vert),
 114 fVertexer(0),
 115 fListOfTracks(0),
 116 fITSclusters(0),
 117 fSixPoints(0),
 118 fOuterStartLayer(0),
 119 fCluLayer(0),
 120 fCluCoord(0)
 121 {
 122   // Standard constructor (Vertex is known and passed to this obj.)
 123   if (geom) {
 124     AliWarning("\"geom\" is actually a dummy argument !");
 125   }
 126   Init();
 127
 128 }
 129
 130 //____________________________________________________________________________
 131 AliITStrackerSA::AliITStrackerSA(const Char_t *geom, AliITSVertexer *vertexer):AliITStrackerMI(0),
 132 fPhiEstimate(0),
 133 fITSStandAlone(0),
 134 fLambdac(0),
 135 fPhic(0),
 136 fCoef1(0),
 137 fCoef2(0),
 138 fCoef3(0),
 139 fNloop(0),
 140 fPhiWin(0),
 141 fLambdaWin(0),
 142 fVert(),
 143 fVertexer(vertexer),
 144 fListOfTracks(0),
 145 fITSclusters(0),
 146 fSixPoints(0),
 147 fOuterStartLayer(0),
 148 fCluLayer(0),
 149 fCluCoord(0)
 150 {
 151   // Standard constructor (Vertex is unknown - vertexer is passed to this obj)
 152   if (geom) {
 153     AliWarning("\"geom\" is actually a dummy argument !");
 154   }
 155   Init();
 156   fVertexer = vertexer;
 157
 158 }
 159
 160 //____________________________________________________________________________
 161 AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA& tracker):AliITStrackerMI(),
 162 fPhiEstimate(tracker.fPhiEstimate),
 163 fITSStandAlone(tracker.fITSStandAlone),
 164 fLambdac(tracker.fLambdac),
 165 fPhic(tracker.fPhic),
 166 fCoef1(tracker.fCoef1),
 167 fCoef2(tracker.fCoef2),
 168 fCoef3(tracker.fCoef3),
 169 fNloop(tracker.fNloop),
 170 fPhiWin(tracker.fPhiWin),
 171 fLambdaWin(tracker.fLambdaWin),
 172 fVert(tracker.fVert),
 173 fVertexer(tracker.fVertexer),
 174 fListOfTracks(tracker.fListOfTracks),
 175 fITSclusters(tracker.fITSclusters),
 176 fSixPoints(tracker.fSixPoints),
 177 fOuterStartLayer(tracker.fOuterStartLayer),
 178 fCluLayer(tracker.fCluLayer),
 179 fCluCoord(tracker.fCluCoord) {
 180   // Copy constructor
 181   for(Int_t i=0;i<2;i++){
 182     fPoint1[i]=tracker.fPoint1[i];
 183     fPoint2[i]=tracker.fPoint2[i];
 184     fPoint3[i]=tracker.fPoint3[i];
 185     fPointc[i]=tracker.fPointc[i];
 186   }
 187   if(tracker.fVertexer && tracker.fVert){
 188     fVert = new AliESDVertex(*tracker.fVert);
 189   }
 190   else {
 191     fVert = tracker.fVert;
 192   }
 193   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 194     fCluLayer[i] = tracker.fCluLayer[i];
 195     fCluCoord[i] = tracker.fCluCoord[i];
 196   }
 197 }
 198 //______________________________________________________________________
 199 AliITStrackerSA& AliITStrackerSA::operator=(const AliITStrackerSA& source){
 200     // Assignment operator.
 201   this->~AliITStrackerSA();
 202   new(this) AliITStrackerSA(source);
 203   return *this;
 204
 205 }
 206
 207 //____________________________________________________________________________
 208 AliITStrackerSA::~AliITStrackerSA(){
 209   // destructor
 210   // if fVertexer is not null, the AliESDVertex obj. is owned by this class
 211   // and is deleted here
 212   if(fVertexer){
 213     if(fVert)delete fVert;
 214   }
 215   fVert = 0;
 216   fVertexer = 0;
 217
 218   if(fPhiWin)delete []fPhiWin;
 219   if(fLambdaWin)delete []fLambdaWin;
 220   fListOfTracks->Delete();
 221   delete fListOfTracks;
 222   if(fCluLayer){
 223     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 224       if(fCluLayer[i]){
 225         fCluLayer[i]->Delete();
 226         delete fCluLayer[i];
 227       }
 228     }
 229     delete [] fCluLayer;
 230   }
 231   if(fCluCoord){
 232     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 233       if(fCluCoord[i]){
 234         fCluCoord[i]->Delete();
 235         delete fCluCoord[i];
 236       }
 237     }
 238     delete [] fCluCoord;
 239   }
 240
 241 }
 242
 243 //____________________________________________________________________________
 244 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
 245 // This method is used to find and fit the tracks. By default the corresponding
 246 // method in the parent class is invoked. In this way a combined tracking
 247 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
 248 // is done in the ITS only. In the standard reconstruction chain this option
 249 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
 250   Int_t rc=0;
 251   if(!fITSStandAlone){
 252     rc=AliITStrackerMI::Clusters2Tracks(event);
 253   }
 254   else {
 255     AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
 256   }
 257   if(!rc) rc=FindTracks(event);
 258   return rc;
 259 }
 260
 261 //____________________________________________________________________________
 262 void AliITStrackerSA::Init(){
 263   //  Reset all data members
 264     fPhiEstimate=0;
 265     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 266     fLambdac=0;
 267     fPhic=0;
 268     fCoef1=0;
 269     fCoef2=0;
 270     fCoef3=0;
 271     fPointc[0]=0;
 272     fPointc[1]=0;
 273     fVert = 0;
 274     fVertexer = 0;
 275     SetWindowSizes();
 276     fITSclusters = 0;
 277     SetSixPoints();
 278     SetOuterStartLayer(0);
 279     SetSAFlag(kFALSE);
 280     fListOfTracks=new TObjArray(0,0);
 281     fCluLayer = 0;
 282     fCluCoord = 0;
 283  }
 284 //_______________________________________________________________________
 285 void AliITStrackerSA::ResetForFinding(){
 286   //  Reset data members used in all loops during track finding
 287     fPhiEstimate=0;
 288     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 289     fLambdac=0;
 290     fPhic=0;
 291     fCoef1=0;
 292     fCoef2=0;
 293     fCoef3=0;
 294     fPointc[0]=0;
 295     fPointc[1]=0;
 296     fListOfTracks->Delete();
 297 }
 298
 299
 300
 301 //______________________________________________________________________
 302 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event){
 303
 304 // Track finder using the ESD object
 305
 306
 307   //controllare numero cluster sui layer1 e 2 (morti?)
 308   //non trova tracce...controllare..
 309
 310   if(!fITSclusters){
 311     Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
 312     return -1;
 313   }
 314
 315
 316   //Reads event and mark clusters of traks already found, with flag kITSin
 317   Int_t nentr=event->GetNumberOfTracks();
 318   while (nentr--) {
 319     AliESDtrack *track=event->GetTrack(nentr);
 320     if (track->GetStatus()&AliESDtrack::kITSin==AliESDtrack::kITSin){
 321       Int_t idx[12];
 322       Int_t ncl = track->GetITSclusters(idx);
 323       for(Int_t k=0;k<ncl;k++){
 324         AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
 325         cll->SetBit(kSAflag);
 326       }
 327     }
 328   }
 329   //Get primary vertex
 330   Double_t primaryVertex[3];
 331   event->GetVertex()->GetXYZ(primaryVertex);
 332   //Creates TClonesArray with clusters for each layer. The clusters already used
 333   //by AliITStrackerMI are not considered
 334   Int_t nclusters[6]={0,0,0,0,0,0};
 335   Int_t dmar[6]={0,0,0,0,0,0};
 336   if (fCluLayer == 0) {
 337     fCluLayer = new TClonesArray*[AliITSgeomTGeo::GetNLayers()];
 338     fCluCoord = new TClonesArray*[AliITSgeomTGeo::GetNLayers()];
 339     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 340       fCluLayer[i]=0;
 341       fCluCoord[i]=0;
 342     }
 343   }
 344   Int_t * firstmod = new Int_t[AliITSgeomTGeo::GetNLayers()];
 345   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 346     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
 347     AliITSlayer &layer=fgLayers[i];
 348     for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 349       AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 350       if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
 351       if(cls->GetQ()==0) continue; //fake clusters dead zones
 352       nclusters[i]++;
 353     }
 354     dmar[i]=0;
 355     delete fCluLayer[i];
 356     fCluLayer[i] = new TClonesArray("AliITSRecPoint",nclusters[i]);
 357     delete fCluCoord[i];
 358     fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
 359   }
 360
 361   for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
 362     TClonesArray &clulay = *fCluLayer[ilay];
 363     TClonesArray &clucoo = *fCluCoord[ilay];
 364     AliITSlayer &layer=fgLayers[ilay];
 365     for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 366       AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 367       if(cls->TestBit(kSAflag)==kTRUE) continue;
 368       if(cls->GetQ()==0) continue;
 369       Double_t phi=0;Double_t lambda=0;
 370       Float_t x=0;Float_t y=0;Float_t z=0;
 371       Float_t sx=0;Float_t sy=0;Float_t sz=0;
 372       GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
 373       GetCoorErrors(cls,sx,sy,sz);
 374       new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
 375       new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
 376       dmar[ilay]++;
 377     }
 378   }
 379
 380   Int_t minNPoints = (fSixPoints ? AliITSgeomTGeo::GetNLayers() : AliITSgeomTGeo::GetNLayers()-1);
 381   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 382     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) {
 383       minNPoints--;
 384     }
 385   }
 386
 387
 388   Int_t ntrack=0;
 389   //loop on the different windows
 390   Int_t * nn = new Int_t[AliITSgeomTGeo::GetNLayers()];//counter for clusters on each layer
 391   for(Int_t nloop=0;nloop<fNloop;nloop++){
 392     for(Int_t ncl=0;ncl<fCluLayer[0]->GetEntries();ncl++){ //loop starting from layer 0
 393
 394       ResetForFinding();
 395       Int_t pflag=0;
 396
 397       AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[0]->At(ncl);
 398
 399       if(!cl) continue;
 400       if (cl->GetQ()<=0) continue;
 401
 402       AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(0,ncl);
 403       fPhic = arr->GetPhi();
 404       fLambdac = arr->GetLambda();
 405       if (TMath::Abs(fLambdac)>0.26*TMath::Pi()) continue;
 406       fPhiEstimate = fPhic;
 407       AliITStrackSA* trs = new AliITStrackSA();
 408       fPoint1[0]=primaryVertex[0];
 409       fPoint1[1]=primaryVertex[1];
 410
 411
 412       fPoint2[0]=arr->GetX();
 413       fPoint2[1]=arr->GetY();
 414       for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){ nn[i]=0;}
 415       nn[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 416       nn[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 417       if(nn[1]>0){
 418         pflag=1;
 419         fPoint3[0] = fPointc[0];
 420         fPoint3[1] = fPointc[1];
 421       }
 422       nn[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 423       if(nn[1]==0 && nn[2]==0) pflag=0;
 424       if(nn[2]!=0 && nn[1]!=0){ pflag=1; UpdatePoints();}
 425       if(nn[2]!=0 && nn[1]==0){
 426         pflag=1;
 427         fPoint3[0]=fPointc[0];
 428         fPoint3[1]=fPointc[1];
 429       }
 430
 431       nn[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 432       pflag=1;
 433       if(nn[3]!=0) UpdatePoints();
 434       nn[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 435       pflag=1;
 436       if(nn[4]!=0) UpdatePoints();
 437       nn[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 438
 439
 440       Int_t layOK=0;
 441       //check of the candidate track
 442       for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers();nnp++) {
 443         if(nn[nnp]!=0) layOK+=1;
 444       }
 445
 446       if(layOK>=minNPoints){
 447         AliITStrackV2* tr2 = 0;
 448         tr2 = FitTrack(trs,primaryVertex);
 449         if(!tr2) continue;
 450
 451         AliESDtrack outtrack;
 452         outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 453         event->AddTrack(&outtrack);
 454         ntrack++;
 455       }
 456       delete trs;
 457      }//end loop on clusters of layer1
 458      //end loop2
 459   }
 460
 461   delete[] nn;
 462
 463
 464   minNPoints--;
 465
 466   //if 5/6 points are required, second loop starting
 467   //from second layer (SPD2), to find tracks with point of
 468   //layer 1 missing
 469   if(!fSixPoints) {
 470     //printf("looking from SPD2\n");
 471     //   counter for clusters on each layer
 472     Int_t * nn = new Int_t[AliITSgeomTGeo::GetNLayers()-1];
 473     for(Int_t nloop=0;nloop<fNloop;nloop++){
 474       Int_t ncl2=fCluLayer[1]->GetEntries();
 475       while(ncl2--){ //loop starting from layer 2
 476         ResetForFinding();
 477         Int_t pflag=0;
 478         AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[1]->At(ncl2);
 479
 480         if(!cl) continue;
 481         AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(1,ncl2);
 482         fPhic = arr->GetPhi();
 483         fLambdac = arr->GetLambda();
 484         fPhiEstimate = fPhic;
 485
 486         AliITStrackSA* trs = new AliITStrackSA();
 487         fPoint1[0]=primaryVertex[0];
 488         fPoint1[1]=primaryVertex[1];
 489
 490         fPoint2[0]=arr->GetX();
 491         fPoint2[1]=arr->GetY();
 492         for(Int_t kk=0;kk<AliITSgeomTGeo::GetNLayers()-1;kk++)nn[kk] = 0;
 493         nn[0] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],
 494                                trs,primaryVertex[2],pflag);
 495         nn[1] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],
 496                                trs,primaryVertex[2],pflag);
 497         if(nn[1]!=0){
 498           pflag=1;
 499           fPoint3[0]=fPointc[0];
 500           fPoint3[1]=fPointc[1];
 501         }
 502         nn[2]= SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],
 503                               trs,primaryVertex[2],pflag);
 504         if(nn[2]!=0){
 505           pflag=1;
 506           UpdatePoints();
 507         }
 508         nn[3]= SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],
 509                               trs,primaryVertex[2],pflag);
 510         if(nn[3]!=0){
 511           pflag=1;
 512           UpdatePoints();
 513         }
 514         nn[4]=SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],
 515                              trs,primaryVertex[2],pflag);
 516
 517         Int_t fl=0;
 518         for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-1;nnp++){
 519           if(nn[nnp]!=0) fl+=1;
 520         }
 521         if(fl>=minNPoints){  // 5/6
 522           AliITStrackV2* tr2 = 0;
 523           tr2 = FitTrack(trs,primaryVertex);
 524           if(!tr2) continue;
 525
 526           AliESDtrack outtrack;
 527           outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 528           event->AddTrack(&outtrack);
 529           ntrack++;
 530
 531         }
 532
 533         delete trs;
 534       }//end loop on clusters of layer2
 535     }
 536
 537     delete [] nn;
 538   } //end opt="5/6"
 539
 540   // search for tracks starting from SPD2, SDD1, SDD2, SSD2
 541   // for cosmics (A. Dainese 31.07.07)
 542   if(fOuterStartLayer>0) {
 543     for(Int_t innLay=1; innLay<=fOuterStartLayer; innLay++) {
 544       //printf("Searching from layer %d outward\n",innLay);
 545       minNPoints=AliITSgeomTGeo::GetNLayers()-innLay;
 546       for(Int_t i=innLay;i<AliITSgeomTGeo::GetNLayers();i++)
 547         if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i))
 548           minNPoints--;
 549       //   counter for clusters on each layer
 550       Int_t * nn = new Int_t[AliITSgeomTGeo::GetNLayers()-innLay];
 551       for(Int_t nloop=0;nloop<fNloop;nloop++){
 552         Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
 553         while(nclInnLay--){ //loop starting from layer innLay
 554           ResetForFinding();
 555           Int_t pflag=0;
 556           AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[innLay]->At(nclInnLay);
 557
 558           if(!cl) continue;
 559           AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(innLay,nclInnLay);
 560           fPhic = arr->GetPhi();
 561           fLambdac = arr->GetLambda();
 562           fPhiEstimate = fPhic;
 563
 564           AliITStrackSA* trs = new AliITStrackSA();
 565           fPoint1[0]=primaryVertex[0];
 566           fPoint1[1]=primaryVertex[1];
 567           fPoint2[0]=arr->GetX();
 568           fPoint2[1]=arr->GetY();
 569
 570           Int_t kk;
 571           for(kk=0;kk<AliITSgeomTGeo::GetNLayers()-innLay;kk++) nn[kk] = 0;
 572
 573           kk=0;
 574           nn[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
 575                                   trs,primaryVertex[2],pflag);
 576           for(Int_t nextLay=innLay+1; nextLay<AliITSgeomTGeo::GetNLayers(); nextLay++) {
 577             kk++;
 578             nn[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 579                                     trs,primaryVertex[2],pflag);
 580             if(nn[kk]!=0){
 581               pflag=1;
 582               if(kk==1) {
 583                 fPoint3[0]=fPointc[0];
 584                 fPoint3[1]=fPointc[1];
 585               } else {
 586                 UpdatePoints();
 587               }
 588             }
 589           }
 590
 591           Int_t fl=0;
 592           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
 593             if(nn[nnp]!=0) fl+=1;
 594           }
 595           if(fl>=minNPoints){
 596             AliITStrackV2* tr2 = 0;
 597             tr2 = FitTrack(trs,primaryVertex);
 598             if(!tr2) continue;
 599
 600
 601             AliESDtrack outtrack;
 602             outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 603             event->AddTrack(&outtrack);
 604             ntrack++;
 605
 606           }
 607
 608           delete trs;
 609         }//end loop on clusters of innLay
 610       } //end loop on window sizes
 611
 612       delete [] nn;
 613     } //end loop on innLay
 614   } //end if(fOuterStartLayer>0)
 615
 616
 617   delete [] firstmod;
 618   Info("FindTracks","Number of found tracks: %d",event->GetNumberOfTracks());
 619   return 0;
 620
 621 }
 622
 623 //________________________________________________________________________
 624
 625 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex) {
 626   //fit of the found track (most general case, also <6 points, layers missing)
 627   // A.Dainese 16.11.07
 628
 629
 630   Int_t * firstmod = new Int_t[AliITSgeomTGeo::GetNLayers()];
 631   TObjArray** listlayer = new TObjArray*[AliITSgeomTGeo::GetNLayers()];
 632   Int_t ** clind=new Int_t*[AliITSgeomTGeo::GetNLayers()];
 633   Int_t ** clmark=new Int_t*[AliITSgeomTGeo::GetNLayers()];
 634   Int_t * nnn = new Int_t[AliITSgeomTGeo::GetNLayers()];
 635   Int_t * kkk = new Int_t[AliITSgeomTGeo::GetNLayers()];
 636
 637   const Int_t maxclu=AliITStrackSA::GetMaxNumberOfClusters();
 638
 639   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 640     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
 641     listlayer[i] = new TObjArray(0,0);
 642     nnn[i]=0;
 643     kkk[i]=0;
 644     clind[i]=new Int_t[maxclu];
 645     clmark[i]=new Int_t[maxclu];
 646     for(Int_t j=0;j<maxclu; j++){
 647       clind[i][j]=0;
 648       clmark[i][j]=0;
 649     }
 650   }
 651
 652
 653   Int_t nclusters = tr->GetNumberOfClustersSA();
 654   for(Int_t ncl=0;ncl<nclusters;ncl++){
 655     Int_t index = tr->GetClusterIndexSA(ncl);
 656     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
 657     if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
 658     Int_t lay = (index & 0xf0000000) >> 28;
 659     listlayer[lay]->AddLast(cl);
 660     Int_t ind=nnn[lay];
 661     clind[lay][ind]=index;
 662   }
 663   delete [] nnn;
 664
 665   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
 666     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
 667       Int_t mark = tr->GetClusterMark(nlay,ncl);
 668       Int_t ind=kkk[nlay];
 669       clmark[nlay][ind]=mark;
 670     }
 671   }
 672
 673   delete [] kkk;
 674
 675   Int_t firstLay=-1,secondLay=-1;
 676   Int_t * end = new Int_t[AliITSgeomTGeo::GetNLayers()];
 677   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 678     if(listlayer[i]->GetEntries()==0) {
 679       end[i]=1;
 680     } else {
 681       end[i]=listlayer[i]->GetEntries();
 682       if(firstLay==-1) {
 683         firstLay=i;
 684       } else if(secondLay==-1) {
 685         secondLay=i;
 686       }
 687     }
 688   }
 689
 690   if(firstLay==-1 || secondLay==-1) {
 691     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 692       delete listlayer[i];
 693       delete clind[i];
 694       delete clmark[i];
 695     }
 696     delete [] listlayer;
 697     delete [] clind;
 698     delete [] clmark;
 699     delete [] firstmod;
 700     delete [] end;
 701     return 0;
 702   }
 703
 704   TClonesArray* listSA = new TClonesArray("AliITStrackSA");
 705   TClonesArray &tri = *listSA;
 706   Int_t nlist=0;
 707
 708   for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
 709     AliITSRecPoint* cl0 = (AliITSRecPoint*)listlayer[0]->At(l0);
 710     for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
 711       AliITSRecPoint* cl1 = (AliITSRecPoint*)listlayer[1]->At(l1);
 712       for(Int_t l2=0;l2<end[2];l2++){  //loop on layer 3
 713         AliITSRecPoint* cl2 = (AliITSRecPoint*)listlayer[2]->At(l2);
 714         for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
 715           AliITSRecPoint* cl3 = (AliITSRecPoint*)listlayer[3]->At(l3);
 716           for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
 717             AliITSRecPoint* cl4 = (AliITSRecPoint*)listlayer[4]->At(l4);
 718             for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
 719               AliITSRecPoint* cl5 = (AliITSRecPoint*)listlayer[5]->At(l5);
 720
 721
 722               Double_t x1,y1,z1,sx1,sy1,sz1;
 723               Double_t x2,y2,z2,sx2,sy2,sz2;
 724               AliITSRecPoint* p1=0;
 725               AliITSRecPoint* p2=0;
 726               Int_t index1=0,index2=0;
 727               Int_t mrk1=0,mrk2=0;
 728
 729               switch(firstLay) {
 730               case 0:
 731                 p1=cl0;
 732                 index1=clind[0][l0];mrk1=clmark[0][l0];
 733                 break;
 734               case 1:
 735                 p1=cl1;
 736                 index1=clind[1][l1];mrk1=clmark[1][l1];
 737                 break;
 738               case 2:
 739                 p1=cl2;
 740                 index1=clind[2][l2];mrk1=clmark[2][l2];
 741                 break;
 742               case 3:
 743                 p1=cl3;
 744                 index1=clind[3][l3];mrk1=clmark[3][l3];
 745                 break;
 746               case 4:
 747                 p1=cl4;
 748                 index1=clind[4][l4];mrk1=clmark[4][l4];
 749                 break;
 750               }
 751
 752               switch(secondLay) {
 753               case 1:
 754                 p2=cl1;
 755                 index2=clind[1][l1];mrk2=clmark[1][l1];
 756                 break;
 757               case 2:
 758                 p2=cl2;
 759                 index2=clind[2][l2];mrk2=clmark[2][l2];
 760                 break;
 761               case 3:
 762                 p2=cl3;
 763                 index2=clind[3][l3];mrk2=clmark[3][l3];
 764                 break;
 765               case 4:
 766                 p2=cl4;
 767                 index2=clind[4][l4];mrk2=clmark[4][l4];
 768                 break;
 769               case 5:
 770                 p2=cl5;
 771                 index2=clind[5][l5];mrk2=clmark[5][l5];
 772                 break;
 773               }
 774
 775               Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
 776
 777               Int_t cln1=mrk1;
 778               Int_t cln2=mrk2;
 779               AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,cln1);
 780               AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,cln2);
 781               x1 = arr1->GetX();
 782               x2 = arr2->GetX();
 783               y1 = arr1->GetY();
 784               y2 = arr2->GetY();
 785               z1 = arr1->GetZ();
 786               z2 = arr2->GetZ();
 787               sx1 = arr1->GetSx();
 788               sx2 = arr2->GetSx();
 789               sy1 = arr1->GetSy();
 790               sy2 = arr2->GetSy();
 791               sz1 = arr1->GetSz();
 792               sz2 = arr2->GetSz();
 793
 794               Int_t layer,ladder,detector;
 795               AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
 796               Float_t yclu1 = p1->GetY();
 797               Float_t zclu1 = p1->GetZ();
 798               Double_t cv=Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
 799               Double_t tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
 800               Double_t phi2 = TMath::ATan2((y2-y1),(x2-x1));
 801               AliITStrackSA* trac = new AliITStrackSA(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
 802
 803
 804               if(cl5!=0) {
 805                 trac->AddClusterV2(5,(clind[5][l5] & 0x0fffffff)>>0);
 806                 trac->AddClusterMark(5,clmark[5][l5]);
 807               }
 808               if(cl4!=0){
 809                 trac->AddClusterV2(4,(clind[4][l4] & 0x0fffffff)>>0);
 810                 trac->AddClusterMark(4,clmark[4][l4]);
 811               }
 812               if(cl3!=0){
 813                 trac->AddClusterV2(3,(clind[3][l3] & 0x0fffffff)>>0);
 814                 trac->AddClusterMark(3,clmark[3][l3]);
 815               }
 816               if(cl2!=0){
 817                 trac->AddClusterV2(2,(clind[2][l2] & 0x0fffffff)>>0);
 818                 trac->AddClusterMark(2,clmark[2][l2]);
 819               }
 820               if(cl1!=0){
 821                 trac->AddClusterV2(1,(clind[1][l1] & 0x0fffffff)>>0);
 822                 trac->AddClusterMark(1,clmark[1][l1]);
 823               }
 824               if(cl0!=0){
 825                 trac->AddClusterV2(0,(clind[0][l0] & 0x0fffffff)>>0);
 826                 trac->AddClusterMark(0,clmark[0][l0]);
 827               }
 828
 829               //fit with Kalman filter using AliITStrackerMI::RefitAt()
 830               AliITStrackMI* ot = new AliITStrackSA(*trac);
 831
 832               ot->ResetCovariance(10.);
 833               ot->ResetClusters();
 834
 835               if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),ot,trac)){ //fit from layer 1 to layer 6
 836                 AliITStrackMI *otrack2 = new AliITStrackMI(*ot);
 837                 otrack2->ResetCovariance(10.);
 838                 otrack2->ResetClusters();
 839                 //fit from layer 6 to layer 1
 840                 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),otrack2,ot)) {
 841                   fListOfTracks->AddLast(otrack2);
 842                   new (tri[nlist]) AliITStrackSA(*trac);
 843                   nlist++;
 844                 } else {
 845                   delete otrack2;
 846                 }
 847
 848               }
 849
 850               delete ot;
 851               delete trac;
 852             }//end loop layer 6
 853           }//end loop layer 5
 854         }//end loop layer 4
 855       }//end loop layer 3
 856     }//end loop layer 2
 857   }//end loop layer 1
 858
 859   delete [] end;
 860
 861
 862
 863   Int_t dim=fListOfTracks->GetEntries();
 864   if(dim==0){
 865     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 866       delete listlayer[i];
 867       delete clind[i];
 868       delete clmark[i];
 869     }
 870     delete [] listlayer;
 871     delete [] clind;
 872     delete [] clmark;
 873     listSA->Delete();
 874     delete listSA;
 875     delete [] firstmod;
 876     return 0;
 877   }
 878
 879   Int_t lowchi2 = FindTrackLowChiSquare(fListOfTracks,dim);
 880   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
 881   AliITStrackSA* trsa = (AliITStrackSA*)listSA->At(lowchi2);
 882
 883   if(otrack==0) {
 884     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 885       delete listlayer[i];
 886       delete clind[i];
 887       delete clmark[i];
 888     }
 889     delete [] listlayer;
 890     delete [] clind;
 891     delete [] clmark;
 892     listSA->Delete();
 893     delete listSA;
 894     delete [] firstmod;
 895     delete [] clind;
 896     delete [] clmark;
 897     return 0;
 898   }
 899   Int_t * indexc = new Int_t[AliITSgeomTGeo::GetNLayers()];
 900   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) indexc[i]=0;
 901   for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
 902     indexc[nind] = otrack->GetClusterIndex(nind);
 903   }
 904   Int_t labl[6][3];
 905   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 906     if(i<otrack->GetNumberOfClusters()) {
 907       AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc[i]);
 908       labl[i][0]=cl->GetLabel(0);
 909       labl[i][1]=cl->GetLabel(1);
 910       labl[i][2]=cl->GetLabel(2);
 911     } else {
 912       labl[i][0]=-1;
 913       labl[i][1]=-1;
 914       labl[i][2]=-1;
 915     }
 916   }
 917   delete [] indexc;
 918
 919   CookLabel(otrack,0.); //MI change - to see fake ratio
 920
 921   Int_t label=FindLabel(labl[0][0],labl[1][0],labl[2][0],labl[3][0],labl[4][0],labl[5][0]);
 922   Int_t lflag=0;
 923   for(Int_t i=0;i<otrack->GetNumberOfClusters();i++)
 924     if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
 925
 926   if(lflag<otrack->GetNumberOfClusters()) label = -label;
 927   otrack->SetLabel(label);
 928
 929   //remove clusters of found track
 930   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
 931     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
 932       Int_t index = trsa->GetClusterMark(nlay,cln);
 933       fCluLayer[nlay]->RemoveAt(index);
 934       RemoveClusterCoord(nlay,index);
 935       fCluLayer[nlay]->Compress();
 936     }
 937   }
 938   listSA->Delete();
 939   delete listSA;
 940
 941   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 942     delete listlayer[i];
 943     delete clind[i];
 944     delete clmark[i];
 945   }
 946   delete [] listlayer;
 947   delete [] firstmod;
 948   delete [] clind;
 949   delete [] clmark;
 950   return otrack;
 951
 952 }
 953
 954
 955
 956 //_______________________________________________________
 957 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
 958   //function used to to find the clusters associated to the track
 959
 960   if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(layer)) return 0;
 961
 962   Int_t nc=0;
 963   AliITSlayer &lay = fgLayers[layer];
 964   Double_t r=lay.GetR();
 965   if(pflag==1){
 966     Float_t cx1,cx2,cy1,cy2;
 967     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
 968     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
 969        return 0;
 970     Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
 971     Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
 972     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
 973   }
 974
 975
 976   Int_t ncl = fCluLayer[layer]->GetEntries();
 977   for (Int_t index=0; index<ncl; index++) {
 978     AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->At(index);
 979     if (!c) continue;
 980     if (c->GetQ()<=0) continue;
 981
 982      AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
 983      Double_t phi = arr->GetPhi();
 984      if (TMath::Abs(phi-fPhiEstimate)>phiwindow) continue;
 985
 986      Double_t lambda = arr->GetLambda();
 987      if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
 988
 989      if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
 990      if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
 991      Int_t orind = arr->GetOrInd();
 992      trs->AddClusterSA(layer,orind);
 993      trs->AddClusterMark(layer,index);
 994
 995      nc++;
 996      fLambdac=lambda;
 997      fPhiEstimate=phi;
 998
 999      fPointc[0]=arr->GetX();
1000      fPointc[1]=arr->GetY();
1001
1002   }
1003   return nc;
1004 }
1005
1006 //________________________________________________________________
1007 void AliITStrackerSA::UpdatePoints(){
1008   //update of points for the estimation of the curvature
1009
1010   fPoint2[0]=fPoint3[0];
1011   fPoint2[1]=fPoint3[1];
1012   fPoint3[0]=fPointc[0];
1013   fPoint3[1]=fPointc[1];
1014
1015
1016 }
1017
1018 //___________________________________________________________________
1019 Int_t AliITStrackerSA::FindEquation(Float_t x1, Float_t y1, Float_t x2, Float_t y2, Float_t x3, Float_t y3,Float_t& a, Float_t& b, Float_t& c){
1020
1021    //given (x,y) of three recpoints (in global coordinates)
1022    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
1023
1024    Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1025    if(den==0) return 0;
1026    a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1027    b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1028    c = -x1*x1-y1*y1-a*x1-b*y1;
1029    return 1;
1030  }
1031 //__________________________________________________________________________
1032  Int_t AliITStrackerSA::FindIntersection(Float_t a1, Float_t b1, Float_t c1, Float_t c2,Float_t& x1,Float_t& y1, Float_t& x2, Float_t& y2){
1033
1034  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
1035  //c2 is -rlayer*rlayer
1036
1037   if(a1==0) return 0;
1038  Double_t m = c2-c1;
1039  Double_t aA = (b1*b1)/(a1*a1)+1;
1040  Double_t bB = (-2*m*b1/(a1*a1));
1041  Double_t cC = c2+(m*m)/(a1*a1);
1042  Double_t dD = bB*bB-4*aA*cC;
1043  if(dD<0) return 0;
1044
1045  y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
1046  y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
1047  x1 = (c2-c1-b1*y1)/a1;
1048  x2 = (c2-c1-b1*y2)/a1;
1049
1050  return 1;
1051 }
1052 //____________________________________________________________________
1053 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
1054 x2,Double_t y2,Double_t x3,Double_t y3){
1055
1056   //calculates the curvature of track
1057   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1058   if(den==0) return 0;
1059   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1060   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1061   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
1062   Double_t xc=-a/2.;
1063
1064   if((a*a+b*b-4*c)<0) return 0;
1065   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
1066   if(rad==0) return 0;
1067
1068   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
1069   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
1070   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
1071   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
1072
1073   return 1/rad;
1074
1075 }
1076
1077
1078 //____________________________________________________________________
1079 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
1080
1081   //Returns the point closest to pp
1082
1083   Double_t diff1 = p1-pp;
1084   Double_t diff2 = p2-pp;
1085
1086   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
1087   else fPhiEstimate=p2;
1088   return fPhiEstimate;
1089
1090 }
1091
1092
1093 //_________________________________________________________________
1094 Int_t AliITStrackerSA::FindTrackLowChiSquare(TObjArray* tracklist, Int_t dim) const {
1095   // returns track with lowers chi square
1096   if(dim==1){
1097     //AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(0);
1098     //return trk;
1099     return 0;
1100   }
1101   if(dim==0) return 0;
1102   Double_t * chi2 = new Double_t[dim];
1103   Int_t * index = new Int_t[dim];
1104   for(Int_t i=0;i<dim;i++){
1105     AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(i);
1106     chi2[i]=trk->GetChi2();
1107     index[i]=i;
1108   }
1109
1110   Int_t w=0;Double_t value;
1111   Int_t lp;
1112   while(w<dim){
1113     for(Int_t j=w+1;j<dim;j++){
1114       if(chi2[w]<chi2[j]){
1115         value=chi2[w];
1116         chi2[w]=chi2[j];
1117         chi2[j]=value;
1118         lp=index[w];
1119         index[w]=index[j];
1120         index[j]=lp;
1121       }
1122     }
1123     w++;
1124   }
1125   Int_t tmp = index[dim-1];
1126   delete [] chi2;
1127   delete [] index;
1128   return tmp;
1129 }
1130
1131 //__________________________________________________________
1132 Int_t AliITStrackerSA::FindLabel(Int_t l0, Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5){
1133
1134   //function used to determine the track label
1135
1136   Int_t lb[6] = {l0,l1,l2,l3,l4,l5};
1137   Int_t aa[6]={1,1,1,1,1,1};
1138   Int_t ff=0;
1139   Int_t ll=0;
1140   Int_t k=0;Int_t w=0;Int_t num=6;
1141   for(Int_t i=5;i>=0;i--) if(lb[i]==-1) num=i;
1142
1143   while(k<num){
1144
1145     for(Int_t i=k+1;i<num;i++){
1146
1147       if(lb[k]==lb[i] && aa[k]!=0){
1148
1149         aa[k]+=1;
1150         aa[i]=0;
1151       }
1152     }
1153     k++;
1154   }
1155
1156   while(w<num){
1157
1158     for(Int_t j=0;j<6;j++){
1159       if(aa[w]<aa[j]){
1160       ff=aa[w];
1161       aa[w]=aa[j];
1162       aa[j]=ff;
1163       ll=lb[w];
1164       lb[w]=lb[j];
1165       lb[j]=ll;
1166      }
1167     }
1168     w++;
1169   }
1170
1171   if(num<1) return -1;
1172   return lb[num-1];
1173 }
1174
1175 //_____________________________________________________________________________
1176 Int_t AliITStrackerSA::Label(Int_t gl1, Int_t gl2, Int_t gl3, Int_t gl4, Int_t gl5, Int_t gl6,Int_t gl7, Int_t gl8, Int_t gl9, Int_t gl10,Int_t gl11,
1177 Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t minNPoints){
1178
1179
1180   //function used to assign label to the found track. If track is fake, the label is negative
1181
1182   Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6};
1183   Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12};
1184   Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18};
1185   Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]);
1186   Int_t lflag=0;Int_t num=6;
1187   if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5;
1188
1189   for(Int_t i=0;i<num;i++){
1190     if(lb0[i]==ll || lb1[i]==ll || lb2[i]==ll) lflag+=1;
1191   }
1192
1193   if(lflag>=minNPoints) return ll;
1194   else return -ll;
1195
1196
1197 }
1198
1199 //_____________________________________________________________________________
1200 void AliITStrackerSA::SetWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1201   // Set sizes of the phi and lambda windows used for track finding
1202   fNloop = n;
1203   if(phi){ // user defined values
1204     fPhiWin = new Double_t[fNloop];
1205     fLambdaWin = new Double_t[fNloop];
1206     for(Int_t k=0;k<fNloop;k++){
1207       fPhiWin[k]=phi[k];
1208       fLambdaWin[k]=lam[k];
1209     }
1210   }
1211   else {  // default values
1212
1213     Double_t phid[33]   = {0.002,0.003,0.004,0.0045,0.0047,
1214                            0.005,0.0053,0.0055,
1215                            0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1216                            0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1217                            0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1218     Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1219                             0.005,0.005,0.006,
1220                             0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1221                             0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1222                             0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1223
1224     if(fNloop!=33){
1225       fNloop = 33;
1226     }
1227
1228
1229     fPhiWin = new Double_t[fNloop];
1230     fLambdaWin = new Double_t[fNloop];
1231
1232     Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
1233
1234     for(Int_t k=0;k<fNloop;k++){
1235       fPhiWin[k]=phid[k]*factor;
1236       fLambdaWin[k]=lambdad[k]*factor;
1237     }
1238
1239   }
1240
1241 }
1242 //_______________________________________________________________________
1243 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1244   //Returns values of phi (azimuthal) and lambda angles for a given cluster
1245 /*
1246   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1247   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1248   Float_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1249
1250   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1251   Double_t phi1=TMath::Pi()/2+alpha;
1252   if (lay==1) phi1+=TMath::Pi();
1253
1254   Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1255   Float_t r=tx*cp+ty*sp;
1256
1257   xyz= r*cp - cl->GetY()*sp;
1258   y= r*sp + cl->GetY()*cp;
1259   z=cl->GetZ();
1260 */
1261   Float_t xyz[3];
1262   cl->GetGlobalXYZ(xyz);
1263   x=xyz[0];
1264   y=xyz[1];
1265   z=xyz[2];
1266
1267   phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1268   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1269 }
1270
1271 //________________________________________________________________________
1272 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1273
1274   //returns sigmax, y, z of cluster in global coordinates
1275 /*
1276   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1277   Int_t lay,lad,det;
1278   AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1279
1280   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1281   Double_t phi=TMath::Pi()/2+alpha;
1282   if (lay==1) phi+=TMath::Pi();
1283
1284   Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1285 */
1286   Float_t covm[6];
1287   cl->GetGlobalCov(covm);
1288   sx=TMath::Sqrt(covm[0]);
1289   sy=TMath::Sqrt(covm[3]);
1290   sz=TMath::Sqrt(covm[5]);
1291 /*
1292   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1293   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1294   sz = TMath::Sqrt(cl->GetSigmaZ2());
1295 */
1296 }