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