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   if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kFALSE) {
 319     while (nentr--) {
 320       AliESDtrack *track=event->GetTrack(nentr);
 321       if (track->GetStatus()&AliESDtrack::kITSin==AliESDtrack::kITSin){
 322         Int_t idx[12];
 323         Int_t ncl = track->GetITSclusters(idx);
 324         for(Int_t k=0;k<ncl;k++){
 325           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
 326           cll->SetBit(kSAflag);
 327         }
 328       }
 329     }
 330   }
 331   //Get primary vertex
 332   Double_t primaryVertex[3];
 333   event->GetVertex()->GetXYZ(primaryVertex);
 334   //Creates TClonesArray with clusters for each layer. The clusters already used
 335   //by AliITStrackerMI are not considered
 336   Int_t nclusters[6]={0,0,0,0,0,0};
 337   Int_t dmar[6]={0,0,0,0,0,0};
 338   if (fCluLayer == 0) {
 339     fCluLayer = new TClonesArray*[AliITSgeomTGeo::GetNLayers()];
 340     fCluCoord = new TClonesArray*[AliITSgeomTGeo::GetNLayers()];
 341     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 342       fCluLayer[i]=0;
 343       fCluCoord[i]=0;
 344     }
 345   }
 346   Int_t * firstmod = new Int_t[AliITSgeomTGeo::GetNLayers()];
 347   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 348     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
 349     AliITSlayer &layer=fgLayers[i];
 350     if (!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) {
 351       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 352         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 353         if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
 354         if(cls->GetQ()==0) continue; //fake clusters dead zones
 355         nclusters[i]++;
 356       }
 357     }
 358     dmar[i]=0;
 359     delete fCluLayer[i];
 360     fCluLayer[i] = new TClonesArray("AliITSRecPoint",nclusters[i]);
 361     delete fCluCoord[i];
 362     fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
 363   }
 364
 365   for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
 366     TClonesArray &clulay = *fCluLayer[ilay];
 367     TClonesArray &clucoo = *fCluCoord[ilay];
 368     AliITSlayer &layer=fgLayers[ilay];
 369     if (!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(ilay)) {
 370       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 371         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 372         if(cls->TestBit(kSAflag)==kTRUE) continue;
 373         if(cls->GetQ()==0) continue;
 374         Double_t phi=0;Double_t lambda=0;
 375         Float_t x=0;Float_t y=0;Float_t z=0;
 376         Float_t sx=0;Float_t sy=0;Float_t sz=0;
 377         GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
 378         GetCoorErrors(cls,sx,sy,sz);
 379         new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
 380         new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
 381         dmar[ilay]++;
 382       }
 383     }
 384   }
 385
 386   Int_t minNPoints = (fSixPoints ? AliITSgeomTGeo::GetNLayers() : AliITSgeomTGeo::GetNLayers()-1);
 387   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 388     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i)) {
 389       minNPoints--;
 390     }
 391   }
 392
 393
 394   Int_t ntrack=0;
 395   //loop on the different windows
 396   Int_t * nn = new Int_t[AliITSgeomTGeo::GetNLayers()];//counter for clusters on each layer
 397   for(Int_t nloop=0;nloop<fNloop;nloop++){
 398     for(Int_t ncl=0;ncl<fCluLayer[0]->GetEntries();ncl++){ //loop starting from layer 0
 399
 400       ResetForFinding();
 401       Int_t pflag=0;
 402
 403       AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[0]->At(ncl);
 404
 405       if(!cl) continue;
 406       if (cl->GetQ()<=0) continue;
 407
 408       AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(0,ncl);
 409       fPhic = arr->GetPhi();
 410       fLambdac = arr->GetLambda();
 411       if (TMath::Abs(fLambdac)>0.26*TMath::Pi()) continue;
 412       fPhiEstimate = fPhic;
 413       AliITStrackSA* trs = new AliITStrackSA();
 414       fPoint1[0]=primaryVertex[0];
 415       fPoint1[1]=primaryVertex[1];
 416
 417
 418       fPoint2[0]=arr->GetX();
 419       fPoint2[1]=arr->GetY();
 420       for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){ nn[i]=0;}
 421       nn[0] = SearchClusters(0,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 422       nn[1] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 423       if(nn[1]>0){
 424         pflag=1;
 425         fPoint3[0] = fPointc[0];
 426         fPoint3[1] = fPointc[1];
 427       }
 428       nn[2] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 429       if(nn[1]==0 && nn[2]==0) pflag=0;
 430       if(nn[2]!=0 && nn[1]!=0){ pflag=1; UpdatePoints();}
 431       if(nn[2]!=0 && nn[1]==0){
 432         pflag=1;
 433         fPoint3[0]=fPointc[0];
 434         fPoint3[1]=fPointc[1];
 435       }
 436
 437       nn[3] = SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 438       pflag=1;
 439       if(nn[3]!=0) UpdatePoints();
 440       nn[4] = SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 441       pflag=1;
 442       if(nn[4]!=0) UpdatePoints();
 443       nn[5] = SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],trs,primaryVertex[2],pflag);
 444
 445
 446       Int_t layOK=0;
 447       //check of the candidate track
 448       for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers();nnp++) {
 449         if(nn[nnp]!=0) layOK+=1;
 450       }
 451
 452       if(layOK>=minNPoints) {
 453         AliITStrackV2* tr2 = 0;
 454         tr2 = FitTrack(trs,primaryVertex);
 455         if(!tr2) continue;
 456
 457         AliESDtrack outtrack;
 458         outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 459         event->AddTrack(&outtrack);
 460         ntrack++;
 461       }
 462       delete trs;
 463      }//end loop on clusters of layer1
 464      //end loop2
 465   }
 466
 467   delete[] nn;
 468
 469
 470   minNPoints--;
 471
 472   //if 5/6 points are required, second loop starting
 473   //from second layer (SPD2), to find tracks with point of
 474   //layer 1 missing
 475   if(!fSixPoints) {
 476     //printf("looking from SPD2\n");
 477     //   counter for clusters on each layer
 478     Int_t * nn = new Int_t[AliITSgeomTGeo::GetNLayers()-1];
 479     for(Int_t nloop=0;nloop<fNloop;nloop++){
 480       Int_t ncl2=fCluLayer[1]->GetEntries();
 481       while(ncl2--){ //loop starting from layer 2
 482         ResetForFinding();
 483         Int_t pflag=0;
 484         AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[1]->At(ncl2);
 485
 486         if(!cl) continue;
 487         AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(1,ncl2);
 488         fPhic = arr->GetPhi();
 489         fLambdac = arr->GetLambda();
 490         fPhiEstimate = fPhic;
 491
 492         AliITStrackSA* trs = new AliITStrackSA();
 493         fPoint1[0]=primaryVertex[0];
 494         fPoint1[1]=primaryVertex[1];
 495
 496         fPoint2[0]=arr->GetX();
 497         fPoint2[1]=arr->GetY();
 498         for(Int_t kk=0;kk<AliITSgeomTGeo::GetNLayers()-1;kk++)nn[kk] = 0;
 499         nn[0] = SearchClusters(1,fPhiWin[nloop],fLambdaWin[nloop],
 500                                trs,primaryVertex[2],pflag);
 501         nn[1] = SearchClusters(2,fPhiWin[nloop],fLambdaWin[nloop],
 502                                trs,primaryVertex[2],pflag);
 503         if(nn[1]!=0){
 504           pflag=1;
 505           fPoint3[0]=fPointc[0];
 506           fPoint3[1]=fPointc[1];
 507         }
 508         nn[2]= SearchClusters(3,fPhiWin[nloop],fLambdaWin[nloop],
 509                               trs,primaryVertex[2],pflag);
 510         if(nn[2]!=0){
 511           pflag=1;
 512           UpdatePoints();
 513         }
 514         nn[3]= SearchClusters(4,fPhiWin[nloop],fLambdaWin[nloop],
 515                               trs,primaryVertex[2],pflag);
 516         if(nn[3]!=0){
 517           pflag=1;
 518           UpdatePoints();
 519         }
 520         nn[4]=SearchClusters(5,fPhiWin[nloop],fLambdaWin[nloop],
 521                              trs,primaryVertex[2],pflag);
 522
 523         Int_t layOK=0;
 524         for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-1;nnp++){
 525           if(nn[nnp]!=0) layOK+=1;
 526         }
 527         if(layOK>=minNPoints){  // 5/6
 528           AliITStrackV2* tr2 = 0;
 529           tr2 = FitTrack(trs,primaryVertex);
 530           if(!tr2) continue;
 531
 532           AliESDtrack outtrack;
 533           outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 534           event->AddTrack(&outtrack);
 535           ntrack++;
 536
 537         }
 538
 539         delete trs;
 540       }//end loop on clusters of layer2
 541     }
 542
 543     delete [] nn;
 544   } //end opt="5/6"
 545
 546   // search for tracks starting from SPD2, SDD1, SDD2, SSD2
 547   // for cosmics (A. Dainese 31.07.07)
 548   if(fOuterStartLayer>0) {
 549     for(Int_t innLay=1; innLay<=fOuterStartLayer; innLay++) {
 550       //printf("Searching from layer %d outward\n",innLay);
 551       minNPoints=AliITSgeomTGeo::GetNLayers()-innLay;
 552       for(Int_t i=innLay;i<AliITSgeomTGeo::GetNLayers();i++)
 553         if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(i))
 554           minNPoints--;
 555       //   counter for clusters on each layer
 556       Int_t * nn = new Int_t[AliITSgeomTGeo::GetNLayers()-innLay];
 557       for(Int_t nloop=0;nloop<fNloop;nloop++){
 558         Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
 559         while(nclInnLay--){ //loop starting from layer innLay
 560           ResetForFinding();
 561           Int_t pflag=0;
 562           AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[innLay]->At(nclInnLay);
 563
 564           if(!cl) continue;
 565           AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(innLay,nclInnLay);
 566           fPhic = arr->GetPhi();
 567           fLambdac = arr->GetLambda();
 568           fPhiEstimate = fPhic;
 569
 570           AliITStrackSA* trs = new AliITStrackSA();
 571           fPoint1[0]=primaryVertex[0];
 572           fPoint1[1]=primaryVertex[1];
 573           fPoint2[0]=arr->GetX();
 574           fPoint2[1]=arr->GetY();
 575
 576           Int_t kk;
 577           for(kk=0;kk<AliITSgeomTGeo::GetNLayers()-innLay;kk++) nn[kk] = 0;
 578
 579           kk=0;
 580           nn[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
 581                                   trs,primaryVertex[2],pflag);
 582           for(Int_t nextLay=innLay+1; nextLay<AliITSgeomTGeo::GetNLayers(); nextLay++) {
 583             kk++;
 584             nn[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 585                                     trs,primaryVertex[2],pflag);
 586             if(nn[kk]!=0){
 587               pflag=1;
 588               if(kk==1) {
 589                 fPoint3[0]=fPointc[0];
 590                 fPoint3[1]=fPointc[1];
 591               } else {
 592                 UpdatePoints();
 593               }
 594             }
 595           }
 596
 597           Int_t layOK=0;
 598           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
 599             if(nn[nnp]!=0) layOK+=1;
 600           }
 601           if(layOK>=minNPoints){
 602             AliITStrackV2* tr2 = 0;
 603             tr2 = FitTrack(trs,primaryVertex);
 604             if(!tr2) continue;
 605
 606
 607             AliESDtrack outtrack;
 608             outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 609             event->AddTrack(&outtrack);
 610             ntrack++;
 611
 612           }
 613
 614           delete trs;
 615         }//end loop on clusters of innLay
 616       } //end loop on window sizes
 617
 618       delete [] nn;
 619     } //end loop on innLay
 620   } //end if(fOuterStartLayer>0)
 621
 622
 623   // search for 1-point tracks, only for cosmics
 624   // (A.Dainese 21.03.08)
 625   if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
 626      TMath::Abs(event->GetMagneticField())<0.01) {
 627     for(Int_t innLay=0; innLay<=fOuterStartLayer; innLay++) {
 628       //   counter for clusters on each layer
 629       Int_t * nn = new Int_t[AliITSgeomTGeo::GetNLayers()-innLay];
 630       for(Int_t nloop=0;nloop<fNloop;nloop++){
 631         Int_t nclInnLay=fCluLayer[innLay]->GetEntries();
 632         while(nclInnLay--){ //loop starting from layer innLay
 633           ResetForFinding();
 634           Int_t pflag=0;
 635           AliITSRecPoint* cl = (AliITSRecPoint*)fCluLayer[innLay]->At(nclInnLay);
 636
 637           if(!cl) continue;
 638           AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(innLay,nclInnLay);
 639           fPhic = arr->GetPhi();
 640           fLambdac = arr->GetLambda();
 641           fPhiEstimate = fPhic;
 642
 643           AliITStrackSA* trs = new AliITStrackSA();
 644           fPoint1[0]=primaryVertex[0];
 645           fPoint1[1]=primaryVertex[1];
 646           fPoint2[0]=arr->GetX();
 647           fPoint2[1]=arr->GetY();
 648
 649           Int_t kk;
 650           for(kk=0;kk<AliITSgeomTGeo::GetNLayers()-innLay;kk++) nn[kk] = 0;
 651
 652           kk=0;
 653           nn[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
 654                                   trs,primaryVertex[2],pflag);
 655           for(Int_t nextLay=innLay+1; nextLay<AliITSgeomTGeo::GetNLayers(); nextLay++) {
 656             kk++;
 657             nn[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 658                                     trs,primaryVertex[2],pflag);
 659             if(nn[kk]!=0){
 660               pflag=1;
 661               if(kk==1) {
 662                 fPoint3[0]=fPointc[0];
 663                 fPoint3[1]=fPointc[1];
 664               } else {
 665                 UpdatePoints();
 666               }
 667             }
 668           }
 669
 670           Int_t layOK=0;
 671           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
 672             if(nn[nnp]!=0) layOK+=1;
 673           }
 674           if(layOK==1) {
 675             AliITStrackV2* tr2 = 0;
 676             Bool_t onePoint = kTRUE;
 677             tr2 = FitTrack(trs,primaryVertex,onePoint);
 678             if(!tr2) continue;
 679
 680             AliESDtrack outtrack;
 681             outtrack.UpdateTrackParams(tr2,AliESDtrack::kITSin);
 682             event->AddTrack(&outtrack);
 683             ntrack++;
 684
 685           }
 686
 687           delete trs;
 688         }//end loop on clusters of innLay
 689       } //end loop on window sizes
 690
 691       delete [] nn;
 692     } //end loop on innLay
 693   } // end search 1-point tracks
 694
 695   delete [] firstmod;
 696   Info("FindTracks","Number of found tracks: %d",event->GetNumberOfTracks());
 697   return 0;
 698
 699 }
 700
 701 //________________________________________________________________________
 702
 703 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
 704   //fit of the found track (most general case, also <6 points, layers missing)
 705   // A.Dainese 16.11.07
 706
 707
 708   Int_t * firstmod = new Int_t[AliITSgeomTGeo::GetNLayers()];
 709   TObjArray** listlayer = new TObjArray*[AliITSgeomTGeo::GetNLayers()];
 710   Int_t ** clind=new Int_t*[AliITSgeomTGeo::GetNLayers()];
 711   Int_t ** clmark=new Int_t*[AliITSgeomTGeo::GetNLayers()];
 712   Int_t * nnn = new Int_t[AliITSgeomTGeo::GetNLayers()];
 713   Int_t * kkk = new Int_t[AliITSgeomTGeo::GetNLayers()];
 714
 715   const Int_t maxclu=AliITStrackSA::GetMaxNumberOfClusters();
 716
 717   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 718     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
 719     listlayer[i] = new TObjArray(0,0);
 720     nnn[i]=0;
 721     kkk[i]=0;
 722     clind[i]=new Int_t[maxclu];
 723     clmark[i]=new Int_t[maxclu];
 724     for(Int_t j=0;j<maxclu; j++){
 725       clind[i][j]=0;
 726       clmark[i][j]=0;
 727     }
 728   }
 729
 730
 731   Int_t nclusters = tr->GetNumberOfClustersSA();
 732   for(Int_t ncl=0;ncl<nclusters;ncl++){
 733     Int_t index = tr->GetClusterIndexSA(ncl);
 734     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
 735     if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
 736     Int_t lay = (index & 0xf0000000) >> 28;
 737     listlayer[lay]->AddLast(cl);
 738     Int_t ind=nnn[lay];
 739     clind[lay][ind]=index;
 740   }
 741   delete [] nnn;
 742
 743   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
 744     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
 745       Int_t mark = tr->GetClusterMark(nlay,ncl);
 746       Int_t ind=kkk[nlay];
 747       clmark[nlay][ind]=mark;
 748     }
 749   }
 750
 751   delete [] kkk;
 752
 753   Int_t firstLay=-1,secondLay=-1;
 754   Int_t * end = new Int_t[AliITSgeomTGeo::GetNLayers()];
 755   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 756     if(listlayer[i]->GetEntries()==0) {
 757       end[i]=1;
 758     } else {
 759       end[i]=listlayer[i]->GetEntries();
 760       if(firstLay==-1) {
 761         firstLay=i;
 762       } else if(secondLay==-1) {
 763         secondLay=i;
 764       }
 765     }
 766   }
 767
 768   if(firstLay==-1 || (secondLay==-1 && !onePoint)) {
 769     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 770       delete listlayer[i];
 771       delete clind[i];
 772       delete clmark[i];
 773     }
 774     delete [] listlayer;
 775     delete [] clind;
 776     delete [] clmark;
 777     delete [] firstmod;
 778     delete [] end;
 779     return 0;
 780   }
 781
 782   TClonesArray* listSA = new TClonesArray("AliITStrackSA");
 783   TClonesArray &tri = *listSA;
 784   Int_t nlist=0;
 785
 786   for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
 787     AliITSRecPoint* cl0 = (AliITSRecPoint*)listlayer[0]->At(l0);
 788     for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
 789       AliITSRecPoint* cl1 = (AliITSRecPoint*)listlayer[1]->At(l1);
 790       for(Int_t l2=0;l2<end[2];l2++){  //loop on layer 3
 791         AliITSRecPoint* cl2 = (AliITSRecPoint*)listlayer[2]->At(l2);
 792         for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
 793           AliITSRecPoint* cl3 = (AliITSRecPoint*)listlayer[3]->At(l3);
 794           for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
 795             AliITSRecPoint* cl4 = (AliITSRecPoint*)listlayer[4]->At(l4);
 796             for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
 797               AliITSRecPoint* cl5 = (AliITSRecPoint*)listlayer[5]->At(l5);
 798
 799
 800               Double_t x1,y1,z1,sx1,sy1,sz1;
 801               Double_t x2,y2,z2,sx2,sy2,sz2;
 802               AliITSRecPoint* p1=0;
 803               AliITSRecPoint* p2=0;
 804               Int_t index1=0,index2=0;
 805               Int_t mrk1=0,mrk2=0;
 806
 807               switch(firstLay) {
 808               case 0:
 809                 p1=cl0;
 810                 index1=clind[0][l0];mrk1=clmark[0][l0];
 811                 break;
 812               case 1:
 813                 p1=cl1;
 814                 index1=clind[1][l1];mrk1=clmark[1][l1];
 815                 break;
 816               case 2:
 817                 p1=cl2;
 818                 index1=clind[2][l2];mrk1=clmark[2][l2];
 819                 break;
 820               case 3:
 821                 p1=cl3;
 822                 index1=clind[3][l3];mrk1=clmark[3][l3];
 823                 break;
 824               case 4:
 825                 p1=cl4;
 826                 index1=clind[4][l4];mrk1=clmark[4][l4];
 827                 break;
 828               }
 829
 830               switch(secondLay) {
 831               case 1:
 832                 p2=cl1;
 833                 index2=clind[1][l1];mrk2=clmark[1][l1];
 834                 break;
 835               case 2:
 836                 p2=cl2;
 837                 index2=clind[2][l2];mrk2=clmark[2][l2];
 838                 break;
 839               case 3:
 840                 p2=cl3;
 841                 index2=clind[3][l3];mrk2=clmark[3][l3];
 842                 break;
 843               case 4:
 844                 p2=cl4;
 845                 index2=clind[4][l4];mrk2=clmark[4][l4];
 846                 break;
 847               case 5:
 848                 p2=cl5;
 849                 index2=clind[5][l5];mrk2=clmark[5][l5];
 850                 break;
 851               default:
 852                 p2=0;
 853                 index2=-1;mrk2=-1;
 854                 break;
 855               }
 856
 857               Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
 858               Int_t layer,ladder,detector;
 859               AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
 860               Float_t yclu1 = p1->GetY();
 861               Float_t zclu1 = p1->GetZ();
 862               Double_t cv=0,tgl2=0,phi2=0;
 863
 864               Int_t cln1=mrk1;
 865               AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,cln1);
 866               x1 = arr1->GetX();
 867               y1 = arr1->GetY();
 868               z1 = arr1->GetZ();
 869               sx1 = arr1->GetSx();
 870               sy1 = arr1->GetSy();
 871               sz1 = arr1->GetSz();
 872
 873               if(secondLay>0) {
 874                 Int_t cln2=mrk2;
 875                 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,cln2);
 876                 x2 = arr2->GetX();
 877                 y2 = arr2->GetY();
 878                 z2 = arr2->GetZ();
 879                 sx2 = arr2->GetSx();
 880                 sy2 = arr2->GetSy();
 881                 sz2 = arr2->GetSz();
 882                 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
 883                 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
 884                 phi2 = TMath::ATan2((y2-y1),(x2-x1));
 885               } else { // special case of 1-point tracks, only for cosmics (B=0)
 886                 x2 = primaryVertex[0];
 887                 y2 = primaryVertex[1];
 888                 z2 = primaryVertex[2];
 889                 cv = 0;
 890                 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 891                 phi2 = TMath::ATan2((y1-y2),(x1-x2));
 892               }
 893
 894
 895               AliITStrackSA* trac = new AliITStrackSA(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
 896
 897
 898               if(cl5!=0) {
 899                 trac->AddClusterV2(5,(clind[5][l5] & 0x0fffffff)>>0);
 900                 trac->AddClusterMark(5,clmark[5][l5]);
 901               }
 902               if(cl4!=0){
 903                 trac->AddClusterV2(4,(clind[4][l4] & 0x0fffffff)>>0);
 904                 trac->AddClusterMark(4,clmark[4][l4]);
 905               }
 906               if(cl3!=0){
 907                 trac->AddClusterV2(3,(clind[3][l3] & 0x0fffffff)>>0);
 908                 trac->AddClusterMark(3,clmark[3][l3]);
 909               }
 910               if(cl2!=0){
 911                 trac->AddClusterV2(2,(clind[2][l2] & 0x0fffffff)>>0);
 912                 trac->AddClusterMark(2,clmark[2][l2]);
 913               }
 914               if(cl1!=0){
 915                 trac->AddClusterV2(1,(clind[1][l1] & 0x0fffffff)>>0);
 916                 trac->AddClusterMark(1,clmark[1][l1]);
 917               }
 918               if(cl0!=0){
 919                 trac->AddClusterV2(0,(clind[0][l0] & 0x0fffffff)>>0);
 920                 trac->AddClusterMark(0,clmark[0][l0]);
 921               }
 922
 923               //fit with Kalman filter using AliITStrackerMI::RefitAt()
 924               AliITStrackMI* ot = new AliITStrackSA(*trac);
 925
 926               ot->ResetCovariance(10.);
 927               ot->ResetClusters();
 928
 929               if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),ot,trac)){ //fit from layer 1 to layer 6
 930                 AliITStrackMI *otrack2 = new AliITStrackMI(*ot);
 931                 otrack2->ResetCovariance(10.);
 932                 otrack2->ResetClusters();
 933                 //fit from layer 6 to layer 1
 934                 if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),otrack2,ot)) {
 935                   fListOfTracks->AddLast(otrack2);
 936                   new (tri[nlist]) AliITStrackSA(*trac);
 937                   nlist++;
 938                 } else {
 939                   delete otrack2;
 940                 }
 941
 942               }
 943
 944               delete ot;
 945               delete trac;
 946             }//end loop layer 6
 947           }//end loop layer 5
 948         }//end loop layer 4
 949       }//end loop layer 3
 950     }//end loop layer 2
 951   }//end loop layer 1
 952
 953   delete [] end;
 954
 955
 956
 957   Int_t dim=fListOfTracks->GetEntries();
 958   if(dim==0){
 959     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 960       delete listlayer[i];
 961       delete clind[i];
 962       delete clmark[i];
 963     }
 964     delete [] listlayer;
 965     delete [] clind;
 966     delete [] clmark;
 967     listSA->Delete();
 968     delete listSA;
 969     delete [] firstmod;
 970     return 0;
 971   }
 972
 973   Int_t lowchi2 = FindTrackLowChiSquare(fListOfTracks,dim);
 974   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
 975   AliITStrackSA* trsa = (AliITStrackSA*)listSA->At(lowchi2);
 976
 977   if(otrack==0) {
 978     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 979       delete listlayer[i];
 980       delete clind[i];
 981       delete clmark[i];
 982     }
 983     delete [] listlayer;
 984     delete [] clind;
 985     delete [] clmark;
 986     listSA->Delete();
 987     delete listSA;
 988     delete [] firstmod;
 989     delete [] clind;
 990     delete [] clmark;
 991     return 0;
 992   }
 993   Int_t * indexc = new Int_t[AliITSgeomTGeo::GetNLayers()];
 994   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) indexc[i]=0;
 995   for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
 996     indexc[nind] = otrack->GetClusterIndex(nind);
 997   }
 998   Int_t labl[6][3];
 999   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
1000     if(i<otrack->GetNumberOfClusters()) {
1001       AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc[i]);
1002       labl[i][0]=cl->GetLabel(0);
1003       labl[i][1]=cl->GetLabel(1);
1004       labl[i][2]=cl->GetLabel(2);
1005     } else {
1006       labl[i][0]=-1;
1007       labl[i][1]=-1;
1008       labl[i][2]=-1;
1009     }
1010   }
1011   delete [] indexc;
1012
1013   CookLabel(otrack,0.); //MI change - to see fake ratio
1014
1015   Int_t label=FindLabel(labl[0][0],labl[1][0],labl[2][0],labl[3][0],labl[4][0],labl[5][0]);
1016   Int_t lflag=0;
1017   for(Int_t i=0;i<otrack->GetNumberOfClusters();i++)
1018     if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
1019
1020   if(lflag<otrack->GetNumberOfClusters()) label = -label;
1021   otrack->SetLabel(label);
1022
1023   //remove clusters of found track
1024   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
1025     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
1026       Int_t index = trsa->GetClusterMark(nlay,cln);
1027       fCluLayer[nlay]->RemoveAt(index);
1028       RemoveClusterCoord(nlay,index);
1029       fCluLayer[nlay]->Compress();
1030     }
1031   }
1032   listSA->Delete();
1033   delete listSA;
1034
1035   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
1036     delete listlayer[i];
1037     delete clind[i];
1038     delete clmark[i];
1039   }
1040   delete [] listlayer;
1041   delete [] firstmod;
1042   delete [] clind;
1043   delete [] clmark;
1044   return otrack;
1045
1046 }
1047
1048
1049
1050 //_______________________________________________________
1051 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
1052   //function used to to find the clusters associated to the track
1053
1054   if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(layer)) return 0;
1055
1056   Int_t nc=0;
1057   AliITSlayer &lay = fgLayers[layer];
1058   Double_t r=lay.GetR();
1059   if(pflag==1){
1060     Float_t cx1,cx2,cy1,cy2;
1061     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
1062     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
1063        return 0;
1064     Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
1065     Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
1066     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
1067   }
1068
1069
1070   Int_t ncl = fCluLayer[layer]->GetEntries();
1071   for (Int_t index=0; index<ncl; index++) {
1072     AliITSRecPoint *c = (AliITSRecPoint*)fCluLayer[layer]->At(index);
1073     if (!c) continue;
1074     if (c->GetQ()<=0) continue;
1075
1076      AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
1077      Double_t phi = arr->GetPhi();
1078      if (TMath::Abs(phi-fPhiEstimate)>phiwindow) continue;
1079
1080      Double_t lambda = arr->GetLambda();
1081      if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
1082
1083      if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
1084      if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
1085      Int_t orind = arr->GetOrInd();
1086      trs->AddClusterSA(layer,orind);
1087      trs->AddClusterMark(layer,index);
1088
1089      nc++;
1090      fLambdac=lambda;
1091      fPhiEstimate=phi;
1092
1093      fPointc[0]=arr->GetX();
1094      fPointc[1]=arr->GetY();
1095
1096   }
1097   return nc;
1098 }
1099
1100 //________________________________________________________________
1101 void AliITStrackerSA::UpdatePoints(){
1102   //update of points for the estimation of the curvature
1103
1104   fPoint2[0]=fPoint3[0];
1105   fPoint2[1]=fPoint3[1];
1106   fPoint3[0]=fPointc[0];
1107   fPoint3[1]=fPointc[1];
1108
1109
1110 }
1111
1112 //___________________________________________________________________
1113 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){
1114
1115    //given (x,y) of three recpoints (in global coordinates)
1116    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
1117
1118    Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1119    if(den==0) return 0;
1120    a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1121    b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1122    c = -x1*x1-y1*y1-a*x1-b*y1;
1123    return 1;
1124  }
1125 //__________________________________________________________________________
1126  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){
1127
1128  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
1129  //c2 is -rlayer*rlayer
1130
1131   if(a1==0) return 0;
1132  Double_t m = c2-c1;
1133  Double_t aA = (b1*b1)/(a1*a1)+1;
1134  Double_t bB = (-2*m*b1/(a1*a1));
1135  Double_t cC = c2+(m*m)/(a1*a1);
1136  Double_t dD = bB*bB-4*aA*cC;
1137  if(dD<0) return 0;
1138
1139  y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
1140  y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
1141  x1 = (c2-c1-b1*y1)/a1;
1142  x2 = (c2-c1-b1*y2)/a1;
1143
1144  return 1;
1145 }
1146 //____________________________________________________________________
1147 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
1148 x2,Double_t y2,Double_t x3,Double_t y3){
1149
1150   //calculates the curvature of track
1151   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
1152   if(den==0) return 0;
1153   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
1154   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
1155   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
1156   Double_t xc=-a/2.;
1157
1158   if((a*a+b*b-4*c)<0) return 0;
1159   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
1160   if(rad==0) return 0;
1161
1162   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
1163   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
1164   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
1165   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
1166
1167   return 1/rad;
1168
1169 }
1170
1171
1172 //____________________________________________________________________
1173 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
1174
1175   //Returns the point closest to pp
1176
1177   Double_t diff1 = p1-pp;
1178   Double_t diff2 = p2-pp;
1179
1180   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
1181   else fPhiEstimate=p2;
1182   return fPhiEstimate;
1183
1184 }
1185
1186
1187 //_________________________________________________________________
1188 Int_t AliITStrackerSA::FindTrackLowChiSquare(TObjArray* tracklist, Int_t dim) const {
1189   // returns track with lowers chi square
1190   if(dim==1){
1191     //AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(0);
1192     //return trk;
1193     return 0;
1194   }
1195   if(dim==0) return 0;
1196   Double_t * chi2 = new Double_t[dim];
1197   Int_t * index = new Int_t[dim];
1198   for(Int_t i=0;i<dim;i++){
1199     AliITStrackV2* trk = (AliITStrackV2*)tracklist->At(i);
1200     chi2[i]=trk->GetChi2();
1201     index[i]=i;
1202   }
1203
1204   Int_t w=0;Double_t value;
1205   Int_t lp;
1206   while(w<dim){
1207     for(Int_t j=w+1;j<dim;j++){
1208       if(chi2[w]<chi2[j]){
1209         value=chi2[w];
1210         chi2[w]=chi2[j];
1211         chi2[j]=value;
1212         lp=index[w];
1213         index[w]=index[j];
1214         index[j]=lp;
1215       }
1216     }
1217     w++;
1218   }
1219   Int_t tmp = index[dim-1];
1220   delete [] chi2;
1221   delete [] index;
1222   return tmp;
1223 }
1224
1225 //__________________________________________________________
1226 Int_t AliITStrackerSA::FindLabel(Int_t l0, Int_t l1, Int_t l2, Int_t l3, Int_t l4, Int_t l5){
1227
1228   //function used to determine the track label
1229
1230   Int_t lb[6] = {l0,l1,l2,l3,l4,l5};
1231   Int_t aa[6]={1,1,1,1,1,1};
1232   Int_t ff=0;
1233   Int_t ll=0;
1234   Int_t k=0;Int_t w=0;Int_t num=6;
1235   for(Int_t i=5;i>=0;i--) if(lb[i]==-1) num=i;
1236
1237   while(k<num){
1238
1239     for(Int_t i=k+1;i<num;i++){
1240
1241       if(lb[k]==lb[i] && aa[k]!=0){
1242
1243         aa[k]+=1;
1244         aa[i]=0;
1245       }
1246     }
1247     k++;
1248   }
1249
1250   while(w<num){
1251
1252     for(Int_t j=0;j<6;j++){
1253       if(aa[w]<aa[j]){
1254       ff=aa[w];
1255       aa[w]=aa[j];
1256       aa[j]=ff;
1257       ll=lb[w];
1258       lb[w]=lb[j];
1259       lb[j]=ll;
1260      }
1261     }
1262     w++;
1263   }
1264
1265   if(num<1) return -1;
1266   return lb[num-1];
1267 }
1268
1269 //_____________________________________________________________________________
1270 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,
1271 Int_t gl12, Int_t gl13, Int_t gl14,Int_t gl15, Int_t gl16, Int_t gl17, Int_t gl18, Int_t minNPoints){
1272
1273
1274   //function used to assign label to the found track. If track is fake, the label is negative
1275
1276   Int_t lb0[6] = {gl1,gl2,gl3,gl4,gl5,gl6};
1277   Int_t lb1[6] = {gl7,gl8,gl9,gl10,gl11,gl12};
1278   Int_t lb2[6] = {gl13,gl14,gl15,gl16,gl17,gl18};
1279   Int_t ll=FindLabel(lb0[0],lb0[1],lb0[2],lb0[3],lb0[4],lb0[5]);
1280   Int_t lflag=0;Int_t num=6;
1281   if(lb0[5]==-1 && lb1[5]==-1 && lb2[5]==-1) num=5;
1282
1283   for(Int_t i=0;i<num;i++){
1284     if(lb0[i]==ll || lb1[i]==ll || lb2[i]==ll) lflag+=1;
1285   }
1286
1287   if(lflag>=minNPoints) return ll;
1288   else return -ll;
1289
1290
1291 }
1292
1293 //_____________________________________________________________________________
1294 void AliITStrackerSA::SetWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1295   // Set sizes of the phi and lambda windows used for track finding
1296   fNloop = n;
1297   if(phi){ // user defined values
1298     fPhiWin = new Double_t[fNloop];
1299     fLambdaWin = new Double_t[fNloop];
1300     for(Int_t k=0;k<fNloop;k++){
1301       fPhiWin[k]=phi[k];
1302       fLambdaWin[k]=lam[k];
1303     }
1304   }
1305   else {  // default values
1306
1307     Double_t phid[33]   = {0.002,0.003,0.004,0.0045,0.0047,
1308                            0.005,0.0053,0.0055,
1309                            0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1310                            0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1311                            0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1312     Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1313                             0.005,0.005,0.006,
1314                             0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1315                             0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1316                             0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1317
1318     if(fNloop!=33){
1319       fNloop = 33;
1320     }
1321
1322
1323     fPhiWin = new Double_t[fNloop];
1324     fLambdaWin = new Double_t[fNloop];
1325
1326     Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
1327
1328     for(Int_t k=0;k<fNloop;k++){
1329       fPhiWin[k]=phid[k]*factor;
1330       fLambdaWin[k]=lambdad[k]*factor;
1331     }
1332
1333   }
1334
1335 }
1336 //_______________________________________________________________________
1337 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Float_t &x, Float_t &y,Float_t &z,Double_t* vertex){
1338   //Returns values of phi (azimuthal) and lambda angles for a given cluster
1339 /*
1340   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1341   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
1342   Float_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);
1343
1344   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1345   Double_t phi1=TMath::Pi()/2+alpha;
1346   if (lay==1) phi1+=TMath::Pi();
1347
1348   Float_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
1349   Float_t r=tx*cp+ty*sp;
1350
1351   xyz= r*cp - cl->GetY()*sp;
1352   y= r*sp + cl->GetY()*cp;
1353   z=cl->GetZ();
1354 */
1355   Float_t xyz[3];
1356   cl->GetGlobalXYZ(xyz);
1357   x=xyz[0];
1358   y=xyz[1];
1359   z=xyz[2];
1360
1361   phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1362   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1363 }
1364
1365 //________________________________________________________________________
1366 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1367
1368   //returns sigmax, y, z of cluster in global coordinates
1369 /*
1370   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1371   Int_t lay,lad,det;
1372   AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1373
1374   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1375   Double_t phi=TMath::Pi()/2+alpha;
1376   if (lay==1) phi+=TMath::Pi();
1377
1378   Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1379 */
1380   Float_t covm[6];
1381   cl->GetGlobalCov(covm);
1382   sx=TMath::Sqrt(covm[0]);
1383   sy=TMath::Sqrt(covm[3]);
1384   sz=TMath::Sqrt(covm[5]);
1385 /*
1386   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1387   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1388   sz = TMath::Sqrt(cl->GetSigmaZ2());
1389 */
1390 }