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 "AliESD.h"
  33 #include "AliESDVertex.h"
  34 #include "AliESDtrack.h"
  35 #include "AliITSVertexer.h"
  36 #include "AliITSclusterTable.h"
  37 #include "AliITSRecPoint.h"
  38 #include "AliITSgeom.h"
  39 #include "AliITStrackSA.h"
  40 #include "AliITStrackerSA.h"
  41 #include "AliRun.h"
  42
  43 ClassImp(AliITStrackerSA)
  44
  45 //____________________________________________________________________________
  46 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(){
  47   // Default constructor
  48   Init();
  49
  50 }
  51 //____________________________________________________________________________
  52 AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom):AliITStrackerMI(geom)
  53 {
  54   // Standard constructor (Vertex is known and passed to this obj.)
  55   Init();
  56   fVert = 0;
  57   fGeom = geom;
  58
  59 }
  60
  61 //____________________________________________________________________________
  62 AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliESDVertex *vert):AliITStrackerMI(geom)
  63 {
  64   // Standard constructor (Vertex is known and passed to this obj.)
  65   Init();
  66   fVert = vert;
  67   fGeom = geom;
  68
  69 }
  70
  71 //____________________________________________________________________________
  72 AliITStrackerSA::AliITStrackerSA(AliITSgeom *geom, AliITSVertexer *vertexer):AliITStrackerMI(geom)
  73 {
  74   // Standard constructor (Vertex is unknown - vertexer is passed to this obj)
  75   Init();
  76   fVertexer = vertexer;
  77   fGeom = geom;
  78
  79 }
  80
  81 //____________________________________________________________________________
  82 AliITStrackerSA::AliITStrackerSA(const AliITStrackerSA& tracker):AliITStrackerMI(){
  83   // Copy constructor
  84   fPhiEstimate = tracker.fPhiEstimate;
  85   for(Int_t i=0;i<2;i++){
  86     fPoint1[i]=tracker.fPoint1[i];
  87     fPoint2[i]=tracker.fPoint2[i];
  88     fPoint3[i]=tracker.fPoint3[i];
  89     fPointc[i]=tracker.fPointc[i];
  90   }
  91   fLambdac = tracker.fLambdac;
  92   fPhic = tracker.fPhic;
  93   fCoef1 = tracker.fCoef1;
  94   fCoef2 = tracker.fCoef2;
  95   fCoef3 = tracker.fCoef3;
  96   fNloop = tracker.fNloop;
  97   fPhiWin = tracker.fPhiWin;
  98   fLambdaWin = tracker.fLambdaWin;
  99   if(tracker.fVertexer && tracker.fVert){
 100     fVert = new AliESDVertex(*tracker.fVert);
 101   }
 102   else {
 103     fVert = tracker.fVert;
 104   }
 105   fVertexer = tracker.fVertexer;
 106   fGeom = tracker.fGeom;
 107   fListOfTracks = tracker.fListOfTracks;
 108   for(Int_t i=0;i<fGeom->GetNlayers();i++){
 109     fCluLayer[i] = tracker.fCluLayer[i];
 110     fCluCoord[i] = tracker.fCluCoord[i];
 111   }
 112 }
 113 //______________________________________________________________________
 114 AliITStrackerSA& AliITStrackerSA::operator=(const AliITStrackerSA& /*source*/){
 115     // Assignment operator. This is a function which is not allowed to be
 116     // done.
 117     Error("operator=","Assignment operator not allowed\n");
 118     return *this;
 119 }
 120
 121 //____________________________________________________________________________
 122 AliITStrackerSA::~AliITStrackerSA(){
 123   // destructor
 124   // if fVertexer is not null, the AliESDVertex obj. is owned by this class
 125   // and is deleted here
 126   if(fVertexer){
 127     if(fVert)delete fVert;
 128   }
 129   fVert = 0;
 130   fVertexer = 0;
 131
 132   if(fPhiWin)delete []fPhiWin;
 133   if(fLambdaWin)delete []fLambdaWin;
 134   fListOfTracks->Delete();
 135   if(fCluLayer){
 136     for(Int_t i=0;i<fGeom->GetNlayers();i++){
 137       if(fCluLayer[i]){
 138         fCluLayer[i]->Delete();
 139         delete fCluLayer[i];
 140       }
 141     }
 142     delete [] fCluLayer;
 143   }
 144   if(fCluCoord){
 145     for(Int_t i=0;i<fGeom->GetNlayers();i++){
 146       if(fCluCoord[i]){
 147         fCluCoord[i]->Delete();
 148         delete fCluCoord[i];
 149       }
 150     }
 151     delete [] fCluCoord;
 152   }
 153
 154   fGeom = 0;
 155 }
 156
 157 //____________________________________________________________________________
 158 void AliITStrackerSA::Init(){
 159   //  Reset all data members
 160     fPhiEstimate=0;
 161     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 162     fLambdac=0;
 163     fPhic=0;
 164     fCoef1=0;
 165     fCoef2=0;
 166     fCoef3=0;
 167     fPointc[0]=0;
 168     fPointc[1]=0;
 169     fVert = 0;
 170     fVertexer = 0;
 171     fGeom = 0;
 172     SetWindowSizes();
 173     fITSclusters = 0;
 174     SetSixPoints();
 175     fListOfTracks=new TObjArray(0,0);
 176     fCluLayer = 0;
 177     fCluCoord = 0;
 178  }
 179 //_______________________________________________________________________
 180 void AliITStrackerSA::ResetForFinding(){
 181   //  Reset data members used in all loops during track finding
 182     fPhiEstimate=0;
 183     for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 184     fLambdac=0;
 185     fPhic=0;
 186     fCoef1=0;
 187     fCoef2=0;
 188     fCoef3=0;
 189     fPointc[0]=0;
 190     fPointc[1]=0;
 191     fListOfTracks->Delete();
 192 }
 193
 194
 195
 196 //______________________________________________________________________
 197 Int_t AliITStrackerSA::FindTracks(AliESD* event){
 198
 199 // Track finder using the ESD object
 200
 201
 202   //controllare numero cluster sui layer1 e 2 (morti?)
 203   //non trova tracce...controllare..
 204
 205   if(!fITSclusters){
 206     Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
 207     return -1;
 208   }
 209
 210
 211   //Reads event and mark clusters of traks already found, with flag kITSin
 212    Int_t nentr=event->GetNumberOfTracks();
 213    while (nentr--) {
 214      AliESDtrack *track=event->GetTrack(nentr);
 215      if (track->GetStatus()&AliESDtrack::kITSin==AliESDtrack::kITSin){
 216        Int_t idx[12];
 217        Int_t ncl = track->GetITSclusters(idx);
 218        for(Int_t k=0;k<ncl;k++){
 219          AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
 220          cll->SetBit(kSAflag);
 221        }
 222      }
 223    }
 224
 225    Double_t primaryVertex[3];
 226    event->GetVertex()->GetXYZ(primaryVertex);
 227    //Creates TClonesArray with clusters for each layer. The clusters already used
 228    //by AliITStrackerMI are not considered
 229
 230    Int_t nclusters[6]={0,0,0,0,0,0};
 231    Int_t dmar[6]={0,0,0,0,0,0};
 232    fCluLayer = new TClonesArray*[fGeom->GetNlayers()];
 233    fCluCoord = new TClonesArray*[fGeom->GetNlayers()];
 234
 235    for(Int_t i=0;i<fGeom->GetNlayers();i++){
 236      AliITSlayer &layer=fgLayers[i];
 237      for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 238        AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 239        if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
 240        if(cls->GetQ()==0) continue; //fake clusters dead zones
 241        nclusters[i]++;
 242      }
 243      dmar[i]=0;
 244      fCluLayer[i] = new TClonesArray("AliITSRecPoint",nclusters[i]);
 245      fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
 246    }
 247
 248
 249    Int_t * firstmod = new Int_t[fGeom->GetNlayers()];
 250    for(Int_t i=0;i<fGeom->GetNlayers();i++){
 251      firstmod[i]=fGeom->GetModuleIndex(i+1,1,1);
 252    }
 253
 254    for(Int_t ilay=0;ilay<fGeom->GetNlayers();ilay++){
 255      TClonesArray &clulay = *fCluLayer[ilay];
 256      TClonesArray &clucoo = *fCluCoord[ilay];
 257      AliITSlayer &layer=fgLayers[ilay];
 258      for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 259        AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 260        if(cls->TestBit(kSAflag)==kTRUE) continue;
 261        if(cls->GetQ()==0) continue;
 262        Double_t phi=0;Double_t lambda=0;
 263        Float_t x=0;Float_t y=0;Float_t z=0;
 264        Float_t sx=0;Float_t sy=0;Float_t sz=0;
 265        Int_t module = cls->GetDetectorIndex()+firstmod[ilay];
 266        GetCoorAngles(cls,module,phi,lambda,x,y,z,primaryVertex);
 267        GetCoorErrors(cls,module,sx,sy,sz);
 268        new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
 269        new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
 270        dmar[ilay]++;
 271      }
 272    }
 273
 274
 275   //Get primary vertex
 276
 277    Int_t ntrack=0;
 278    //loop on the different windows
 279    for(Int_t nloop=0;nloop<fNloop;nloop++){
 280      for(Int_t ncl=0;ncl<fCluLayer[0]->GetEntries();ncl++){ //loop starting from layer 0
 281
 282        ResetForFinding();
 283        Int_t pflag=0;
 284
 285        AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[0]->At(ncl);
 286
 287        if(!cl) continue;
 288
 289
 290        if (cl->GetQ()<=0) continue;
 291
 292        AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(0,ncl);
 293        fPhic = arr->GetPhi();
 294        fLambdac = arr->GetLambda();
 295        if (TMath::Abs(fLambdac)>0.26*TMath::Pi()) continue;
 296        fPhiEstimate = fPhic;
 297        AliITStrackSA* trs = new AliITStrackSA();
 298        fPoint1[0]=primaryVertex[0];
 299        fPoint1[1]=primaryVertex[1];
 300
 301
 302        fPoint2[0]=arr->GetX();
 303        fPoint2[1]=arr->GetY();
 304        Int_t * nn = new Int_t[fGeom->GetNlayers()];//counter for clusters on each layer
 305        for(Int_t i=0;i<fGeom->GetNlayers();i++){ nn[i]=0;}
 306        nn[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 307
 308        nn[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 309       if(nn[1]>0){
 310         pflag=1;
 311         fPoint3[0] = fPointc[0];
 312         fPoint3[1] = fPointc[1];
 313       }
 314       nn[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 315       if(nn[1]==0 && nn[2]==0) pflag=0;
 316       if(nn[2]!=0 && nn[1]!=0){ pflag=1; UpdatePoints();}
 317       if(nn[2]!=0 && nn[1]==0){
 318         pflag=1;
 319         fPoint3[0]=fPointc[0];
 320         fPoint3[1]=fPointc[1];
 321       }
 322
 323       nn[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 324       pflag=1;
 325       if(nn[3]!=0) UpdatePoints();
 326       nn[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 327       pflag=1;
 328       if(nn[4]!=0) UpdatePoints();
 329       nn[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 330
 331
 332       Int_t layOK=0;
 333       Int_t numberofpoints;
 334       if(fSixPoints) numberofpoints=6;  //check of the candidate track
 335       else numberofpoints=5;           //if track is good (with the required number
 336       for(Int_t nnp=0;nnp<fGeom->GetNlayers();nnp++){    //of points) it is written on file
 337         if(nn[nnp]!=0) layOK+=1;
 338       }
 339       if(layOK>=numberofpoints){
 340
 341         AliITStrackV2* tr2 = FitTrack(trs,primaryVertex);
 342
 343         if(tr2==0) continue;
 344
 345         AliESDtrack outtrack;
 346         outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 347         event->AddTrack(&outtrack);
 348         ntrack++;
 349
 350       }
 351
 352
 353       delete trs;
 354       delete[] nn;
 355
 356      }//end loop on clusters of layer1
 357
 358      //end loop2
 359    }
 360
 361    //if 5/6 points are required, second loop starting
 362    //from second layer, to find tracks with point of
 363    //layer 1 missing
 364
 365    if(!fSixPoints){
 366      //   counter for clusters on each layer
 367      Int_t * nn = new Int_t[fGeom->GetNlayers()-1];
 368      for(Int_t nloop=0;nloop<fNloop;nloop++){
 369        Int_t ncl2=fCluLayer[1]->GetEntries();
 370        while(ncl2--){ //loop starting from layer 2
 371          ResetForFinding();
 372          Int_t pflag=0;
 373          AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[1]->At(ncl2);
 374
 375          if(!cl) continue;
 376          AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(1,ncl2);
 377          fPhic = arr->GetPhi();
 378          fLambdac = arr->GetLambda();
 379          fPhiEstimate = fPhic;
 380
 381          AliITStrackSA* trs = new AliITStrackSA();
 382          fPoint1[0]=primaryVertex[0];
 383          fPoint1[1]=primaryVertex[1];
 384
 385          fPoint2[0]=arr->GetX();
 386          fPoint2[1]=arr->GetY();
 387          for(Int_t kk=0;kk<fGeom->GetNlayers()-1;kk++)nn[kk] = 0;
 388          nn[0] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],
 389                                    trs,primaryVertex[2],pflag);
 390          nn[1] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],
 391                                    trs,primaryVertex[2],pflag);
 392          if(nn[1]!=0){
 393            pflag=1;
 394            fPoint3[0]=fPointc[0];
 395            fPoint3[1]=fPointc[1];
 396          }
 397          nn[2]= SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],
 398                                    trs,primaryVertex[2],pflag);
 399          if(nn[2]!=0){
 400            pflag=1;
 401            UpdatePoints();
 402          }
 403          nn[3]= SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],
 404                                    trs,primaryVertex[2],pflag);
 405          if(nn[3]!=0){
 406            pflag=1;
 407            UpdatePoints();
 408          }
 409          nn[4]=SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],
 410                                    trs,primaryVertex[2],pflag);
 411
 412          Int_t fl=0;
 413          for(Int_t nnp=0;nnp<fGeom->GetNlayers()-1;nnp++){
 414            if(nn[nnp]!=0) fl+=1;
 415          }
 416          if(fl>=5){  // 5/6
 417            AliITStrackV2* tr2 = FitTrack(trs,primaryVertex);
 418            if(tr2==0){
 419              continue;
 420            }
 421
 422            AliESDtrack outtrack;
 423            outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 424            event->AddTrack(&outtrack);
 425            ntrack++;
 426
 427          }
 428
 429
 430         delete trs;
 431       }//end loop on clusters of layer2
 432      }
 433
 434     delete [] nn;
 435   }  //end opt="5/6"
 436
 437   delete [] firstmod;
 438   Info("FindTracks","Number of found tracks: %d",event->GetNumberOfTracks());
 439   return 0;
 440
 441 }
 442
 443
 444
 445
 446
 447 //________________________________________________________________________
 448
 449 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex){
 450   //fit of the found track
 451
 452
 453   Int_t * firstmod = new Int_t[fGeom->GetNlayers()];
 454   for(Int_t i=0;i<fGeom->GetNlayers();i++){
 455     firstmod[i]=fGeom->GetModuleIndex(i+1,1,1);
 456   }
 457
 458   Int_t nclusters = tr->GetNumberOfClustersSA();
 459   TObjArray** listlayer = new TObjArray*[fGeom->GetNlayers()];
 460   for(Int_t i=0;i<fGeom->GetNlayers();i++){
 461     listlayer[i] = new TObjArray(0,0);
 462   }
 463
 464   TArrayI clind0(20);
 465   TArrayI clind1(20);
 466   TArrayI clind2(20);
 467   TArrayI clind3(20);
 468   TArrayI clind4(20);
 469   TArrayI clind5(20);
 470
 471   TArrayI mark0(20);
 472   TArrayI mark1(20);
 473   TArrayI mark2(20);
 474   TArrayI mark3(20);
 475   TArrayI mark4(20);
 476   TArrayI mark5(20);
 477
 478
 479   Int_t * nnn = new Int_t[fGeom->GetNlayers()];
 480   Int_t * kkk = new Int_t[fGeom->GetNlayers()];
 481   for(Int_t i=0;i<fGeom->GetNlayers();i++) {nnn[i]=0;kkk[i]=0;}
 482
 483   for(Int_t ncl=0;ncl<nclusters;ncl++){
 484     Int_t index = tr->GetClusterIndexSA(ncl);
 485     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
 486     if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
 487     Int_t lay = (index & 0xf0000000) >> 28;
 488     if(lay==0) { listlayer[0]->AddLast(cl); clind0[nnn[0]]=index;nnn[0]++;}
 489     if(lay==1) { listlayer[1]->AddLast(cl); clind1[nnn[1]]=index;nnn[1]++;}
 490     if(lay==2) { listlayer[2]->AddLast(cl); clind2[nnn[2]]=index;nnn[2]++;}
 491     if(lay==3) { listlayer[3]->AddLast(cl); clind3[nnn[3]]=index;nnn[3]++;}
 492     if(lay==4) { listlayer[4]->AddLast(cl); clind4[nnn[4]]=index;nnn[4]++;}
 493     if(lay==5) { listlayer[5]->AddLast(cl); clind5[nnn[5]]=index;nnn[5]++;}
 494   }
 495   delete [] nnn;
 496
 497   for(Int_t nlay=0;nlay<fGeom->GetNlayers();nlay++){
 498     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
 499       Int_t mark = tr->GetClusterMark(nlay,ncl);
 500       if(nlay==0) { mark0[kkk[0]]=mark;kkk[0]++;}
 501       if(nlay==1) { mark1[kkk[1]]=mark;kkk[1]++;}
 502       if(nlay==2) { mark2[kkk[2]]=mark;kkk[2]++;}
 503       if(nlay==3) { mark3[kkk[3]]=mark;kkk[3]++;}
 504       if(nlay==4) { mark4[kkk[4]]=mark;kkk[4]++;}
 505       if(nlay==5) { mark5[kkk[5]]=mark;kkk[5]++;}
 506
 507     }
 508   }
 509
 510   delete [] kkk;
 511
 512
 513   Int_t * end = new Int_t[fGeom->GetNlayers()];
 514   for(Int_t i=0;i<fGeom->GetNlayers();i++){
 515     if(listlayer[i]->GetEntries()==0) end[i]=1;
 516     else end[i]=listlayer[i]->GetEntries();
 517   }
 518
 519   TClonesArray* listSA = new TClonesArray("AliITStrackSA");
 520   TClonesArray &tri = *listSA;
 521   Int_t nlist=0;
 522
 523   if(end[0]==0) end[0]=1; //for tracks with cluster on layer 0 missing
 524   for(Int_t l1=0;l1<end[0];l1++){//loop on layer 1
 525     AliITSRecPoint* cl0 = (AliITSRecPoint*)listlayer[0]->At(l1);
 526     Double_t x1,y1,z1,sx1,sy1,sz1;
 527     Double_t x2,y2,z2,sx2,sy2,sz2;
 528     AliITSRecPoint* p1=0;
 529     AliITSRecPoint* p2=0;
 530     Int_t index1=clind0[l1];
 531     Int_t index2=0;
 532     Int_t mrk1 = mark0[l1];
 533     Int_t mrk2 = 0;
 534     Int_t lay1=0;
 535     Int_t lay2=1;
 536     Int_t module1=-1;
 537     for(Int_t l2=0;l2<end[1];l2++){//loop on layer 2
 538       AliITSRecPoint* cl1 = (AliITSRecPoint*)listlayer[1]->At(l2);
 539       index2=clind1[l2];
 540       mrk2 = mark1[l2];
 541       for(Int_t l3=0;l3<end[2];l3++){  //loop on layer 3
 542         AliITSRecPoint* cl2 = (AliITSRecPoint*)listlayer[2]->At(l3);
 543
 544         if(cl0==0 && cl1!=0) {
 545           p2 = cl2;index1=clind2[l3];mrk1=mark2[l3];lay1=2;
 546           p1=cl1;
 547           module1 = p1->GetDetectorIndex()+firstmod[1];
 548         }
 549         if(cl0!=0 && cl1==0){
 550           p1=cl0;
 551           p2=cl2;index2=clind2[l3];mrk2=mark2[l3];lay2=2;
 552           module1 = p1->GetDetectorIndex()+firstmod[0];
 553         }
 554         if(cl0!=0 && cl1!=0){
 555           p1=cl0;
 556           p2=cl1;
 557           module1 = p1->GetDetectorIndex()+firstmod[0];
 558         }
 559
 560         Int_t cln1=mrk1;
 561         Int_t cln2=mrk2;
 562         AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay1,cln1);
 563         AliITSclusterTable* arr1= (AliITSclusterTable*)GetClusterCoord(lay2,cln2);
 564         x1 = arr->GetX();
 565         x2 = arr1->GetX();
 566         y1 = arr->GetY();
 567         y2 = arr1->GetY();
 568         z1 = arr->GetZ();
 569         z2 = arr1->GetZ();
 570         sx1 = arr->GetSx();
 571         sx2 = arr1->GetSx();
 572         sy1 = arr->GetSy();
 573         sy2 = arr1->GetSy();
 574         sz1 = arr->GetSz();
 575         sz2 = arr1->GetSz();
 576
 577         Int_t layer,ladder,detector;
 578         fGeom->GetModuleId(module1,layer,ladder,detector);
 579         Float_t yclu1 = p1->GetY();
 580         Float_t zclu1 = p1->GetZ();
 581         Double_t cv=Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
 582
 583         Double_t tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
 584         Double_t phi2 = TMath::ATan2((y2-y1),(x2-x1));
 585
 586         for(Int_t l4=0;l4<end[3];l4++){ //loop on layer 4
 587           AliITSRecPoint* cl3 = (AliITSRecPoint*)listlayer[3]->At(l4);
 588           for(Int_t l5=0;l5<end[4];l5++){ //loop on layer 5
 589             AliITSRecPoint* cl4 = (AliITSRecPoint*)listlayer[4]->At(l5);
 590             for(Int_t l6=0;l6<end[5];l6++){ //loop on layer 6
 591               AliITSRecPoint* cl5 = (AliITSRecPoint*)listlayer[5]->At(l6);
 592               AliITStrackSA* trac = new AliITStrackSA(fGeom,layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
 593
 594               if(cl5!=0) {
 595                 trac->AddClusterV2(5,(clind5[l6] & 0x0fffffff)>>0);
 596                 trac->AddClusterMark(5,mark5[l6]);
 597               }
 598               if(cl4!=0){
 599                 trac->AddClusterV2(4,(clind4[l5] & 0x0fffffff)>>0);
 600                 trac->AddClusterMark(4,mark4[l5]);
 601               }
 602               if(cl3!=0){
 603                 trac->AddClusterV2(3,(clind3[l4] & 0x0fffffff)>>0);
 604                 trac->AddClusterMark(3,mark3[l4]);
 605               }
 606               if(cl2!=0){
 607                 trac->AddClusterV2(2,(clind2[l3] & 0x0fffffff)>>0);
 608                 trac->AddClusterMark(2,mark2[l3]);
 609               }
 610               if(cl1!=0){
 611                 trac->AddClusterV2(1,(clind1[l2] & 0x0fffffff)>>0);
 612                 trac->AddClusterMark(1,mark1[l2]);
 613               }
 614               if(cl0!=0){
 615                 trac->AddClusterV2(0,(clind0[l1] & 0x0fffffff)>>0);
 616                 trac->AddClusterMark(0,mark0[l1]);
 617               }
 618               //fit with Kalman filter using AliITStrackerMI::RefitAt()
 619
 620
 621               AliITStrackMI* ot = new AliITStrackSA(*trac);
 622
 623               ot->ResetCovariance();
 624               ot->ResetClusters();
 625
 626               if(RefitAt(49.,ot,trac)){ //fit from layer 1 to layer 6
 627                 AliITStrackMI *otrack2 = new AliITStrackMI(*ot);
 628                 otrack2->ResetCovariance();
 629                 otrack2->ResetClusters();
 630                 //fit from layer 6 to layer 1
 631                 if(RefitAt(3.7,otrack2,ot)) {
 632                   fListOfTracks->AddLast(otrack2);
 633                   new (tri[nlist]) AliITStrackSA(*trac);
 634                   nlist++;
 635                 } else {
 636                   delete otrack2;
 637                 }
 638
 639               }
 640
 641               delete ot;
 642               delete trac;
 643             }//end loop layer 6
 644           }//end loop layer 5
 645         }//end loop layer 4
 646
 647       }//end loop layer 3
 648     }//end loop layer 2
 649   }//end loop layer 1
 650
 651   delete [] end;
 652
 653
 654
 655   Int_t dim=fListOfTracks->GetEntries();
 656   if(dim==0){
 657     for(Int_t i=0;i<fGeom->GetNlayers();i++){
 658       delete listlayer[i];
 659     }
 660     delete [] listlayer;
 661     listSA->Delete();
 662     delete listSA;
 663     delete [] firstmod;
 664     return 0;
 665   }
 666
 667   Int_t lowchi2 = FindTrackLowChiSquare(fListOfTracks,dim);
 668   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
 669   AliITStrackSA* trsa = (AliITStrackSA*)listSA->At(lowchi2);
 670
 671   if(otrack==0) {
 672     for(Int_t i=0;i<fGeom->GetNlayers();i++){
 673       delete listlayer[i];
 674     }
 675     delete [] listlayer;
 676     listSA->Delete();
 677     delete listSA;
 678     delete [] firstmod;
 679     return 0;
 680   }
 681   Int_t * indexc = new Int_t[fGeom->GetNlayers()];
 682   for(Int_t i=0;i<fGeom->GetNlayers();i++) indexc[i]=0;
 683   for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
 684     indexc[nind] = otrack->GetClusterIndex(nind);
 685   }
 686   AliITSRecPoint* cl0 = (AliITSRecPoint*)GetCluster(indexc[0]);
 687   AliITSRecPoint* cl1 = (AliITSRecPoint*)GetCluster(indexc[1]);
 688   AliITSRecPoint* cl2 = (AliITSRecPoint*)GetCluster(indexc[2]);
 689   AliITSRecPoint* cl3 = (AliITSRecPoint*)GetCluster(indexc[3]);
 690   AliITSRecPoint* cl4 = (AliITSRecPoint*)GetCluster(indexc[4]);
 691   Int_t labl[3]={-1,-1,-1};
 692   if(otrack->GetNumberOfClusters()==fGeom->GetNlayers()){
 693     AliITSRecPoint* cl5 = (AliITSRecPoint*)GetCluster(indexc[5]);
 694     labl[0]=cl5->GetLabel(0);
 695     labl[1]=cl5->GetLabel(1);
 696     labl[2]=cl5->GetLabel(2);
 697   }
 698   delete [] indexc;
 699   if(otrack->GetNumberOfClusters()==(fGeom->GetNlayers()-1)){
 700     labl[0]=-1;
 701     labl[1]=-1;
 702     labl[2]=-1;
 703   }
 704   Int_t numberofpoints;
 705   if(fSixPoints) numberofpoints=6;
 706   else numberofpoints=5;
 707   CookLabel(otrack,0.); //MI change - to see fake ratio
 708   Int_t label =  Label(cl0->GetLabel(0),cl1->GetLabel(0),
 709                        cl2->GetLabel(0),cl3->GetLabel(0),
 710                        cl4->GetLabel(0),labl[0],
 711                        cl0->GetLabel(1),cl1->GetLabel(1),
 712                        cl2->GetLabel(1),cl3->GetLabel(1),
 713                        cl4->GetLabel(1),labl[1],
 714                        cl0->GetLabel(2),cl1->GetLabel(2),
 715                        cl2->GetLabel(2),cl3->GetLabel(2),
 716                        cl4->GetLabel(2),labl[2],numberofpoints);
 717
 718   otrack->SetLabel(label);
 719   //remove clusters of found track
 720
 721   for(Int_t nlay=0;nlay<fGeom->GetNlayers();nlay++){
 722     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
 723       Int_t index = trsa->GetClusterMark(nlay,cln);
 724       fCluLayer[nlay]->RemoveAt(index);
 725       RemoveClusterCoord(nlay,index);
 726       fCluLayer[nlay]->Compress();
 727     }
 728   }
 729   listSA->Delete();
 730   delete listSA;
 731
 732   for(Int_t i=0;i<fGeom->GetNlayers();i++){
 733     delete listlayer[i];
 734   }
 735   delete [] listlayer;
 736   delete [] firstmod;
 737   return otrack;
 738
 739 }
 740
 741
 742
 743 //_______________________________________________________
 744 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
 745   //function used to to find the clusters associated to the track
 746   Int_t nc=0;
 747   AliITSlayer &lay = fgLayers[layer];
 748   Double_t r=lay.GetR();
 749    if(pflag==1){
 750     Float_t cx1,cx2,cy1,cy2;
 751     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
 752     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
 753        return 0;
 754     Double_t fi1=TMath::ATan2(cy1,cx1);
 755     Double_t fi2=TMath::ATan2(cy2,cx2);
 756     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
 757   }
 758
 759
 760   Int_t ncl = fCluLayer[layer]->GetEntries();
 761   for (Int_t index=0; index<ncl; index++) {
 762     AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->At(index);
 763     if (!c) continue;
 764     if (c->GetQ()<=0) continue;
 765
 766      AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
 767      Double_t phi = arr->GetPhi();
 768      if (TMath::Abs(phi-fPhiEstimate)>phiwindow) continue;
 769
 770      Double_t lambda = arr->GetLambda();
 771      if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
 772
 773      if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
 774      if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
 775      Int_t orind = arr->GetOrInd();
 776      trs->AddClusterSA(layer,orind);
 777      trs->AddClusterMark(layer,index);
 778
 779      nc++;
 780      fLambdac=lambda;
 781      fPhiEstimate=phi;
 782
 783      fPointc[0]=arr->GetX();
 784      fPointc[1]=arr->GetY();
 785
 786   }
 787   return nc;
 788 }
 789
 790 //________________________________________________________________
 791 void AliITStrackerSA::UpdatePoints(){
 792   //update of points for the estimation of the curvature
 793
 794   fPoint2[0]=fPoint3[0];
 795   fPoint2[1]=fPoint3[1];
 796   fPoint3[0]=fPointc[0];
 797   fPoint3[1]=fPointc[1];
 798
 799
 800 }
 801
 802 //___________________________________________________________________
 803 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){
 804
 805    //given (x,y) of three recpoints (in global coordinates)
 806    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
 807
 808    Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 809    if(den==0) return 0;
 810    a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 811    b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 812    c = -x1*x1-y1*y1-a*x1-b*y1;
 813    return 1;
 814  }
 815 //__________________________________________________________________________
 816  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){
 817
 818  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
 819  //c2 is -rlayer*rlayer
 820
 821   if(a1==0) return 0;
 822  Double_t m = c2-c1;
 823  Double_t aA = (b1*b1)/(a1*a1)+1;
 824  Double_t bB = (-2*m*b1/(a1*a1));
 825  Double_t cC = c2+(m*m)/(a1*a1);
 826  Double_t dD = bB*bB-4*aA*cC;
 827  if(dD<0) return 0;
 828
 829  y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
 830  y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
 831  x1 = (c2-c1-b1*y1)/a1;
 832  x2 = (c2-c1-b1*y2)/a1;
 833
 834  return 1;
 835 }
 836 //____________________________________________________________________
 837 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
 838 x2,Double_t y2,Double_t x3,Double_t y3){
 839
 840   //calculates the curvature of track
 841   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 842   if(den==0) return 0;
 843   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 844   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 845   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
 846   Double_t xc=-a/2.;
 847
 848   if((a*a+b*b-4*c)<0) return 0;
 849   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
 850   if(rad==0) return 0;
 851
 852   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
 853   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
 854   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
 855   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
 856
 857   return 1/rad;
 858
 859 }
 860
 861
 862 //____________________________________________________________________
 863 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
 864
 865   //Returns the point closest to pp
 866
 867   Double_t diff1 = p1-pp;
 868   Double_t diff2 = p2-pp;
 869
 870   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
 871   else fPhiEstimate=p2;
 872   return fPhiEstimate;
 873
 874 }
 875
 876
 877 //_________________________________________________________________
 878 Int_t AliITStrackerSA::FindTrackLowChiSquare(TObjArray* tracklist, Int_t dim) const {
 879   // returns track with lowes chi square
 880   if(dim==1){
 881     //AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(0);
 882     //return trk;
 883     return 0;
 884   }
 885   //if(dim==0) return 0;
 886   Double_t * chi2 = new Double_t[dim];
 887   Int_t * index = new Int_t[dim];
 888   for(Int_t i=0;i<dim;i++){
 889     AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(i);
 890     chi2[i]=trk->GetChi2();
 891     index[i]=i;
 892   }
 893
 894   Int_t w=0;Double_t value;
 895   Int_t lp;
 896   while(w<dim){
 897     for(Int_t j=w+1;j<dim;j++){
 898       if(chi2[w]<chi2[j]){
 899         value=chi2[w];
 900         chi2[w]=chi2[j];
 901         chi2[j]=value;
 902         lp=index[w];
 903         index[w]=index[j];
 904         index[j]=lp;
 905       }
 906     }
 907     w++;
 908   }
 909
 910   delete [] chi2;
 911   delete [] index;
 912   return index[dim-1];
 913 }
 914
 915 //__________________________________________________________
 916 Int_t AliITStrackerSA::FindLabel(Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5, Int_t l6){
 917
 918   //function used to determine the track label
 919
 920   Int_t lb[6] = {l1,l2,l3,l4,l5,l6};
 921   Int_t aa[6]={1,1,1,1,1,1};
 922   Int_t ff=0;
 923   Int_t ll=0;
 924   Int_t k=0;Int_t w=0;Int_t num=6;
 925   if(lb[5]==-1) num=5;
 926
 927   while(k<num){
 928
 929     for(Int_t i=k+1;i<num;i++){
 930
 931       if(lb[k]==lb[i] && aa[k]!=0){
 932
 933         aa[k]+=1;
 934         aa[i]=0;
 935       }
 936     }
 937   k++;
 938   }
 939
 940   while(w<num){
 941
 942     for(Int_t j=0;j<6;j++){
 943       if(aa[w]<aa[j]){
 944       ff=aa[w];
 945       aa[w]=aa[j];
 946       aa[j]=ff;
 947       ll=lb[w];
 948       lb[w]=lb[j];
 949       lb[j]=ll;
 950      }
 951     }
 952   w++;
 953   }
 954   if(num==6)  return lb[5];
 955   else return lb[4];
 956 }
 957
 958 //_____________________________________________________________________________
 959 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,
 960 Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t numberofpoints){
 961
 962
 963   //function used to assign label to the found track. If track is fake, the label is negative
 964
 965   Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6};
 966   Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12};
 967   Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18};
 968   Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]);
 969   Int_t lflag=0;Int_t num=6;
 970   if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5;
 971
 972   for(Int_t i=0;i<num;i++){
 973     if(lb0[i]==ll || lb1[i]==ll || lb2[i]==ll) lflag+=1;
 974   }
 975
 976   if(lflag>=numberofpoints) return ll;
 977   else return -ll;
 978
 979
 980 }
 981
 982 //_____________________________________________________________________________
 983 void AliITStrackerSA::SetWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
 984   // Set sizes of the phi and lambda windows used for track finding
 985   fNloop = n;
 986   if(phi){ // user defined values
 987     fPhiWin = new Double_t[fNloop];
 988     fLambdaWin = new Double_t[fNloop];
 989     for(Int_t k=0;k<fNloop;k++){
 990       fPhiWin[k]=phi[k];
 991       fLambdaWin[k]=lam[k];
 992     }
 993   }
 994   else {  // default values
 995
 996     Double_t phid[33]   = {0.002,0.003,0.004,0.0045,0.0047,
 997                            0.005,0.0053,0.0055,
 998                            0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
 999                            0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1000                            0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1001     Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1002                             0.005,0.005,0.006,
1003                             0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1004                             0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1005                             0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1006
1007     if(fNloop!=33){
1008       fNloop = 33;
1009     }
1010
1011
1012     fPhiWin = new Double_t[fNloop];
1013     fLambdaWin = new Double_t[fNloop];
1014
1015     for(Int_t k=0;k<fNloop;k++){
1016       fPhiWin[k]=phid[k];
1017       fLambdaWin[k]=lambdad[k];
1018     }
1019
1020   }
1021
1022 }
1023 //_______________________________________________________________________
1024 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Int_t module,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1025   //Returns values of phi (azimuthal) and lambda angles for a given cluster
1026
1027   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1028   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1029   Float_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1030
1031   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1032   Double_t phi1=TMath::Pi()/2+alpha;
1033   if (lay==1) phi1+=TMath::Pi();
1034
1035   Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1036   Float_t r=tx*cp+ty*sp;
1037
1038   x= r*cp - cl->GetY()*sp;
1039   y= r*sp + cl->GetY()*cp;
1040   z=cl->GetZ();
1041
1042   phi=TMath::ATan2(y,x);
1043   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1044 }
1045
1046 //________________________________________________________________________
1047 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl, Int_t module,Float_t &sx,Float_t &sy, Float_t &sz){
1048
1049   //returns x,y,z of cluster in global coordinates
1050
1051   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1052   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1053
1054   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1055   Double_t phi=TMath::Pi()/2+alpha;
1056   if (lay==1) phi+=TMath::Pi();
1057
1058   Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1059
1060   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1061   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1062   sz = TMath::Sqrt(cl->GetSigmaZ2());
1063
1064 }