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