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