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