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