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