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