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