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 ITS tracker class                        //
  20 //  Origin:  Elisabetta Crescio - crescio@to.infn.it     //
  21 //  Updated: Francesco Prino    - prino@to.infn.it       //
  22 ///////////////////////////////////////////////////////////
  23
  24 #include <stdlib.h>
  25
  26 #include <TArrayI.h>
  27 #include <TBranch.h>
  28 #include <TObjArray.h>
  29 #include <TTree.h>
  30
  31 #include "AliESDEvent.h"
  32 #include "AliESDVertex.h"
  33 #include "AliESDtrack.h"
  34 #include "AliITSVertexer.h"
  35 #include "AliITSclusterTable.h"
  36 #include "AliITSRecPoint.h"
  37 #include "AliITSgeomTGeo.h"
  38 #include "AliITStrackSA.h"
  39 #include "AliITStrackerSA.h"
  40 #include "AliITSReconstructor.h"
  41 #include "AliLog.h"
  42 #include "AliRun.h"
  43
  44 ClassImp(AliITStrackerSA)
  45
  46 //____________________________________________________________________________
  47 AliITStrackerSA::AliITStrackerSA():AliITStrackerMI(),
  48 fPhiEstimate(0),
  49 fITSStandAlone(0),
  50 fLambdac(0),
  51 fPhic(0),
  52 fCoef1(0),
  53 fCoef2(0),
  54 fCoef3(0),
  55 fNloop(0),
  56 fPhiWin(0),
  57 fLambdaWin(0),
  58 fListOfTracks(0),
  59 fListOfSATracks(0),
  60 fITSclusters(0),
  61 fInwardFlag(0),
  62 fOuterStartLayer(0),
  63 fInnerStartLayer(5),
  64 fMinNPoints(0),
  65 fMinQ(0.),
  66 fCluCoord(0){
  67   // Default constructor
  68   Init();
  69
  70 }
  71 //____________________________________________________________________________
  72 AliITStrackerSA::AliITStrackerSA(const Char_t *geom):AliITStrackerMI(0),
  73 fPhiEstimate(0),
  74 fITSStandAlone(0),
  75 fLambdac(0),
  76 fPhic(0),
  77 fCoef1(0),
  78 fCoef2(0),
  79 fCoef3(0),
  80 fNloop(0),
  81 fPhiWin(0),
  82 fLambdaWin(0),
  83 fListOfTracks(0),
  84 fListOfSATracks(0),
  85 fITSclusters(0),
  86 fInwardFlag(0),
  87 fOuterStartLayer(0),
  88 fInnerStartLayer(5),
  89 fMinNPoints(0),
  90 fMinQ(0.),
  91 fCluCoord(0)
  92 {
  93   // Standard constructor (Vertex is known and passed to this obj.)
  94   if (geom) {
  95     AliWarning("\"geom\" is actually a dummy argument !");
  96   }
  97
  98   Init();
  99
 100 }
 101
 102 //____________________________________________________________________________
 103 AliITStrackerSA::~AliITStrackerSA(){
 104   // destructor
 105
 106   if(fPhiWin)delete []fPhiWin;
 107   if(fLambdaWin)delete []fLambdaWin;
 108   fListOfTracks->Delete();
 109   delete fListOfTracks;
 110   fListOfSATracks->Delete();
 111   delete fListOfSATracks;
 112   if(fCluCoord){
 113     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 114       if(fCluCoord[i]){
 115         fCluCoord[i]->Delete();
 116         delete fCluCoord[i];
 117       }
 118     }
 119     delete [] fCluCoord;
 120   }
 121 }
 122
 123 //____________________________________________________________________________
 124 Int_t AliITStrackerSA::Clusters2Tracks(AliESDEvent *event){
 125 // This method is used to find and fit the tracks. By default the corresponding
 126 // method in the parent class is invoked. In this way a combined tracking
 127 // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
 128 // is done in the ITS only. In the standard reconstruction chain this option
 129 // can be set via AliReconstruction::SetOption("ITS","onlyITS")
 130   Int_t rc=0;
 131
 132   if(!fITSStandAlone){
 133     rc=AliITStrackerMI::Clusters2Tracks(event);
 134   }
 135   else {
 136     AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
 137   }
 138   if(!rc){
 139     rc=FindTracks(event,kFALSE);
 140     Int_t nSPDcontr=0;
 141     const AliESDVertex *spdv = event->GetPrimaryVertexSPD();
 142     if(spdv) nSPDcontr = spdv->GetNContributors();
 143     if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE &&
 144        nSPDcontr<=AliITSReconstructor::GetRecoParam()->GetMaxSPDcontrForSAToUseAllClusters()) {
 145       rc=FindTracks(event,kTRUE);
 146     }
 147   }
 148   return rc;
 149 }
 150
 151 //____________________________________________________________________________
 152 void AliITStrackerSA::Init(){
 153   //  Reset all data members
 154     fPhiEstimate=0;
 155     for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 156     fLambdac=0;
 157     fPhic=0;
 158     fCoef1=0;
 159     fCoef2=0;
 160     fCoef3=0;
 161     fPointc[0]=0;
 162     fPointc[1]=0;
 163     Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
 164     if(nLoops==33){
 165       SetFixedWindowSizes();
 166     }else{
 167       Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
 168       Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
 169       Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
 170       Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
 171       SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
 172     }
 173     fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
 174     fITSclusters = 0;
 175     SetOuterStartLayer(1);
 176     SetSAFlag(kFALSE);
 177     fListOfTracks=new TClonesArray("AliITStrackMI",100);
 178     fListOfSATracks=new TClonesArray("AliITStrackSA",100);
 179     fCluCoord = 0;
 180     fMinNPoints = 3;
 181  }
 182 //_______________________________________________________________________
 183 void AliITStrackerSA::ResetForFinding(){
 184   //  Reset data members used in all loops during track finding
 185     fPhiEstimate=0;
 186     for(Int_t i=0;i<2;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 187     fLambdac=0;
 188     fPhic=0;
 189     fCoef1=0;
 190     fCoef2=0;
 191     fCoef3=0;
 192     fPointc[0]=0;
 193     fPointc[1]=0;
 194     fListOfTracks->Clear();
 195     fListOfSATracks->Clear();
 196 }
 197
 198
 199
 200 //______________________________________________________________________
 201 Int_t AliITStrackerSA::FindTracks(AliESDEvent* event, Bool_t useAllClusters){
 202
 203 // Track finder using the ESD object
 204
 205   AliDebug(2,Form(" field is %f",event->GetMagneticField()));
 206   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
 207
 208   if(!fITSclusters){
 209     Fatal("FindTracks","ITS cluster tree is not accessed!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
 210     return -1;
 211   }
 212   //Reads event and mark clusters of traks already found, with flag kITSin
 213   Int_t nentr=event->GetNumberOfTracks();
 214   if(!useAllClusters) {
 215     while (nentr--) {
 216       AliESDtrack *track=event->GetTrack(nentr);
 217       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
 218         Int_t idx[12];
 219         Int_t ncl = track->GetITSclusters(idx);
 220         for(Int_t k=0;k<ncl;k++){
 221           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
 222           cll->SetBit(kSAflag);
 223         }
 224       }
 225     }
 226   }else{
 227     while (nentr--) {
 228       AliESDtrack *track=event->GetTrack(nentr);
 229       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
 230         Int_t idx[12];
 231         Int_t ncl = track->GetITSclusters(idx);
 232         for(Int_t k=0;k<ncl;k++){
 233           AliITSRecPoint* cll = (AliITSRecPoint*)GetCluster(idx[k]);
 234           cll->ResetBit(kSAflag);
 235         }
 236       }
 237     }
 238   }
 239   //Get primary vertex
 240   Double_t primaryVertex[3];
 241   event->GetVertex()->GetXYZ(primaryVertex);
 242   //Creates TClonesArray with clusters for each layer. The clusters already used
 243   //by AliITStrackerMI are not considered
 244   Int_t nclusters[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
 245   Int_t dmar[AliITSgeomTGeo::kNLayers]={0,0,0,0,0,0};
 246   if (fCluCoord == 0) {
 247     fCluCoord = new TClonesArray*[AliITSgeomTGeo::kNLayers];
 248     for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 249       fCluCoord[i]=0;
 250     }
 251   }
 252   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++){
 253     AliITSlayer &layer=fgLayers[i];
 254     if (!ForceSkippingOfLayer(i)) {
 255       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 256         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 257         if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
 258         if(cls->GetQ()==0) continue; //fake clusters dead zones
 259         if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
 260         nclusters[i]++;
 261       }
 262     }
 263     dmar[i]=0;
 264     if(!fCluCoord[i]){
 265       fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
 266     }else{
 267       fCluCoord[i]->Delete();
 268       fCluCoord[i]->Expand(nclusters[i]);
 269     }
 270   }
 271
 272   for(Int_t ilay=0;ilay<AliITSgeomTGeo::GetNLayers();ilay++){
 273     TClonesArray &clucoo = *fCluCoord[ilay];
 274     AliITSlayer &layer=fgLayers[ilay];
 275     if (!ForceSkippingOfLayer(ilay)) {
 276       for(Int_t cli=0;cli<layer.GetNumberOfClusters();cli++){
 277         AliITSRecPoint* cls = (AliITSRecPoint*)layer.GetCluster(cli);
 278         if(cls->TestBit(kSAflag)==kTRUE) continue;
 279         if(cls->GetQ()==0) continue;
 280         if(ilay>1 && cls->GetQ()<=fMinQ) continue;
 281         Double_t phi=0;Double_t lambda=0;
 282         Double_t x=0;Double_t y=0;Double_t z=0;
 283         Double_t sx=0;Double_t sy=0;Double_t sz=0;
 284         GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
 285         GetCoorErrors(cls,sx,sy,sz);
 286         new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
 287         dmar[ilay]++;
 288       }
 289     }
 290   }
 291
 292   // track counter
 293   Int_t ntrack=0;
 294
 295   static Int_t nClusLay[AliITSgeomTGeo::kNLayers];//counter for clusters on each layer
 296   Int_t startLayForSeed=0;
 297   Int_t lastLayForSeed=fOuterStartLayer;
 298   Int_t nSeedSteps=lastLayForSeed-startLayForSeed;
 299   Int_t seedStep=1;
 300   if(fInwardFlag){
 301     startLayForSeed=AliITSgeomTGeo::GetNLayers()-1;
 302     lastLayForSeed=fInnerStartLayer;
 303     nSeedSteps=startLayForSeed-lastLayForSeed;
 304     seedStep=-1;
 305   }
 306
 307   // loop on minimum number of points
 308   for(Int_t iMinNPoints=AliITSgeomTGeo::GetNLayers(); iMinNPoints>=fMinNPoints; iMinNPoints--) {
 309
 310     // loop on starting layer for track finding
 311     for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) {
 312       Int_t theLay=startLayForSeed+iSeedLay*seedStep;
 313       if(ForceSkippingOfLayer(theLay)) continue;
 314       Int_t minNPoints=iMinNPoints-theLay;
 315       if(fInwardFlag) minNPoints=iMinNPoints-(AliITSgeomTGeo::GetNLayers()-1-theLay);
 316       for(Int_t i=theLay+1;i<AliITSgeomTGeo::GetNLayers();i++)
 317         if(ForceSkippingOfLayer(i))
 318           minNPoints--;
 319       if(minNPoints<fMinNPoints) continue;
 320
 321       // loop on phi and lambda window size
 322       for(Int_t nloop=0;nloop<fNloop;nloop++){
 323         Int_t nclTheLay=fCluCoord[theLay]->GetEntries();
 324         while(nclTheLay--){
 325           ResetForFinding();
 326           Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex);
 327           if(!useRP) continue;
 328           AliITStrackSA trs;
 329
 330           Int_t pflag=0;
 331           Int_t kk;
 332           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
 333
 334           kk=0;
 335           nClusLay[kk] = SearchClusters(theLay,fPhiWin[nloop],fLambdaWin[nloop],
 336                                         &trs,primaryVertex[2],pflag);
 337           Int_t nextLay=theLay+seedStep;
 338           Bool_t goon=kTRUE;
 339           if(nextLay<0 || nextLay == 6) goon = kFALSE;
 340           while(goon){
 341             kk++;
 342             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 343                                             &trs,primaryVertex[2],pflag);
 344             if(nClusLay[kk]!=0){
 345               pflag=1;
 346               if(kk==1) {
 347                 fPoint3[0]=fPointc[0];
 348                 fPoint3[1]=fPointc[1];
 349               } else {
 350                 UpdatePoints();
 351               }
 352             }
 353             nextLay+=seedStep;
 354             if(nextLay<0 || nextLay==6) goon=kFALSE;
 355           }
 356
 357
 358           Int_t layOK=0;
 359           if(!fInwardFlag){
 360             for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-theLay;nnp++){
 361               if(nClusLay[nnp]!=0) layOK+=1;
 362             }
 363           }else{
 364             for(Int_t nnp=theLay; nnp>=0; nnp--){
 365               if(nClusLay[nnp]!=0) layOK+=1;
 366             }
 367           }
 368           if(layOK>=minNPoints){
 369             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]));
 370             AliITStrackV2* tr2 = 0;
 371             tr2 = FitTrack(&trs,primaryVertex);
 372             if(!tr2){
 373               continue;
 374             }
 375             AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
 376
 377             StoreTrack(tr2,event,useAllClusters);
 378             ntrack++;
 379
 380           }
 381
 382         }//end loop on clusters of theLay
 383       } //end loop on window sizes
 384     } //end loop on theLay
 385   }//end loop on min points
 386
 387   // search for 1-point tracks in SPD, only for cosmics
 388   // (A.Dainese 21.03.08)
 389   if(AliITSReconstructor::GetRecoParam()->GetSAOnePointTracks() &&
 390      TMath::Abs(event->GetMagneticField())<0.01) {
 391     Int_t outerLayer=1; // only SPD
 392     for(Int_t innLay=0; innLay<=TMath::Min(1,fOuterStartLayer); innLay++) {
 393       //   counter for clusters on each layer
 394
 395       for(Int_t nloop=0;nloop<fNloop;nloop++){
 396         Int_t nclInnLay=fCluCoord[innLay]->GetEntries();
 397         while(nclInnLay--){ //loop starting from layer innLay
 398           ResetForFinding();
 399           Bool_t useRP=SetFirstPoint(innLay,nclInnLay,primaryVertex);
 400           if(!useRP) continue;
 401           AliITStrackSA trs;
 402
 403           Int_t pflag=0;
 404           Int_t kk;
 405           for(kk=0;kk<AliITSgeomTGeo::GetNLayers();kk++) nClusLay[kk] = 0;
 406
 407           kk=0;
 408           nClusLay[kk] = SearchClusters(innLay,fPhiWin[nloop],fLambdaWin[nloop],
 409                                   &trs,primaryVertex[2],pflag);
 410           for(Int_t nextLay=innLay+1; nextLay<=outerLayer; nextLay++) {
 411             kk++;
 412             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 413                                     &trs,primaryVertex[2],pflag);
 414             if(nClusLay[kk]!=0){
 415               pflag=1;
 416               if(kk==1) {
 417                 fPoint3[0]=fPointc[0];
 418                 fPoint3[1]=fPointc[1];
 419               } else {
 420                 UpdatePoints();
 421               }
 422             }
 423           }
 424
 425           Int_t layOK=0;
 426           for(Int_t nnp=0;nnp<AliITSgeomTGeo::GetNLayers()-innLay;nnp++){
 427             if(nClusLay[nnp]!=0) layOK+=1;
 428           }
 429           if(layOK==1) {
 430             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]));
 431             AliITStrackV2* tr2 = 0;
 432             Bool_t onePoint = kTRUE;
 433             tr2 = FitTrack(&trs,primaryVertex,onePoint);
 434             if(!tr2){
 435               continue;
 436             }
 437             AliDebug(2,Form("----NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
 438
 439             StoreTrack(tr2,event,useAllClusters);
 440             ntrack++;
 441
 442           }
 443
 444         }//end loop on clusters of innLay
 445       } //end loop on window sizes
 446
 447     } //end loop on innLay
 448   } // end search 1-point tracks
 449
 450   if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
 451   ResetForFinding();
 452   return 0;
 453
 454 }
 455
 456 //________________________________________________________________________
 457
 458 AliITStrackV2* AliITStrackerSA::FitTrack(AliITStrackSA* tr,Double_t *primaryVertex,Bool_t onePoint) {
 459   //fit of the found track (most general case, also <6 points, layers missing)
 460   // A.Dainese 16.11.07
 461
 462
 463   const Int_t kMaxClu=AliITStrackSA::kMaxNumberOfClusters;
 464
 465   static Int_t firstmod[AliITSgeomTGeo::kNLayers];
 466   static Int_t clind[AliITSgeomTGeo::kNLayers][kMaxClu];
 467   static Int_t clmark[AliITSgeomTGeo::kNLayers][kMaxClu];
 468   static Int_t end[AliITSgeomTGeo::kNLayers];
 469   static Int_t indices[AliITSgeomTGeo::kNLayers];
 470
 471   static AliITSRecPoint *listlayer[AliITSgeomTGeo::kNLayers][kMaxClu];
 472
 473   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 474     firstmod[i]=AliITSgeomTGeo::GetModuleIndex(i+1,1,1);
 475     end[i]=0;
 476     for(Int_t j=0;j<kMaxClu; j++){
 477       clind[i][j]=0;
 478       clmark[i][j]=0;
 479       listlayer[i][j]=0;
 480    }
 481   }
 482
 483
 484   Int_t nclusters = tr->GetNumberOfClustersSA();
 485   for(Int_t ncl=0;ncl<nclusters;ncl++){
 486     Int_t index = tr->GetClusterIndexSA(ncl);
 487     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(index);
 488     Int_t lay = (index & 0xf0000000) >> 28;
 489     Int_t nInLay=end[lay];
 490     listlayer[lay][nInLay]=cl;
 491     clind[lay][nInLay]=index;
 492     end[lay]++;
 493   }
 494
 495   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
 496     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
 497       Int_t mark = tr->GetClusterMark(nlay,ncl);
 498       clmark[nlay][ncl]=mark;
 499     }
 500   }
 501
 502
 503   Int_t firstLay=-1,secondLay=-1;
 504   for(Int_t i=0;i<AliITSgeomTGeo::GetNLayers();i++) {
 505     if(end[i]==0) {
 506       end[i]=1;
 507     }else{
 508       if(firstLay==-1) {
 509         firstLay=i;
 510       } else if(secondLay==-1) {
 511         secondLay=i;
 512       }
 513     }
 514   }
 515
 516   if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
 517   TClonesArray &arrMI= *fListOfTracks;
 518   TClonesArray &arrSA= *fListOfSATracks;
 519   Int_t nFoundTracks=0;
 520
 521
 522   for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
 523     indices[0]=l0;
 524     for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
 525       indices[1]=l1;
 526       for(Int_t l2=0;l2<end[2];l2++){  //loop on layer 3
 527         indices[2]=l2;
 528         for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
 529           indices[3]=l3;
 530           for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
 531             indices[4]=l4;
 532             for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
 533               indices[5]=l5;
 534
 535               // estimate curvature from 2 innermost points (or innermost point + vertex)
 536
 537               Int_t iFirstLay=indices[firstLay];
 538               Int_t mrk1=clmark[firstLay][iFirstLay];
 539
 540               AliITSRecPoint* p1=(AliITSRecPoint*)listlayer[firstLay][iFirstLay];
 541               Int_t module1 = p1->GetDetectorIndex()+firstmod[firstLay];
 542               Int_t layer,ladder,detector;
 543               AliITSgeomTGeo::GetModuleId(module1,layer,ladder,detector);
 544               Double_t yclu1 = p1->GetY();
 545               Double_t zclu1 = p1->GetZ();
 546
 547               Double_t x1,y1,z1;
 548               Double_t x2,y2,z2;
 549               Double_t cv=0,tgl2=0,phi2=0;
 550               AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,mrk1);
 551               x1 = arr1->GetX();
 552               y1 = arr1->GetY();
 553               z1 = arr1->GetZ();
 554
 555               if(secondLay>0) {
 556                 Int_t iSecondLay=indices[secondLay];
 557                 Int_t mrk2=clmark[secondLay][iSecondLay];
 558                 AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,mrk2);
 559                 x2 = arr2->GetX();
 560                 y2 = arr2->GetY();
 561                 z2 = arr2->GetZ();
 562                 cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
 563                 tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
 564                 phi2 = TMath::ATan2((y2-y1),(x2-x1));
 565               } else { // special case of 1-point tracks, only for cosmics (B=0)
 566                 x2 = primaryVertex[0];
 567                 y2 = primaryVertex[1];
 568                 z2 = primaryVertex[2];
 569                 cv = 0;
 570                 tgl2 = (z1-z2)/TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 571                 phi2 = TMath::ATan2((y1-y2),(x1-x2));
 572               }
 573
 574               // create track and attach it the RecPoints
 575               AliITStrackSA trac(layer,ladder,detector,yclu1,zclu1,phi2,tgl2,cv,1);
 576               for(Int_t iLay=5; iLay>=0; iLay--){
 577                 Int_t iInLay=indices[iLay];
 578                 AliITSRecPoint* cl=(AliITSRecPoint*)listlayer[iLay][iInLay];
 579                 if(cl!=0){
 580                   trac.AddClusterV2(iLay,(clind[iLay][iInLay] & 0x0fffffff)>>0);
 581                   trac.AddClusterMark(iLay,clmark[iLay][iInLay]);
 582                 }
 583               }
 584
 585               //fit with Kalman filter using AliITStrackerMI::RefitAt()
 586               AliITStrackSA ot(trac);
 587
 588               ot.ResetCovariance(10.);
 589               ot.ResetClusters();
 590
 591               // Propagate inside the innermost layer with a cluster
 592               if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
 593
 594                 if(RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&ot,&trac)){ //fit from layer 1 to layer 6
 595                   AliITStrackMI otrack2(ot);
 596                   otrack2.ResetCovariance(10.);
 597                   otrack2.ResetClusters();
 598                   //fit from layer 6 to layer 1
 599                   if(RefitAt(AliITSRecoParam::GetrInsideSPD1(),&otrack2,&ot)) {
 600                     new(arrMI[nFoundTracks]) AliITStrackMI(otrack2);
 601                     new(arrSA[nFoundTracks]) AliITStrackSA(trac);
 602                     ++nFoundTracks;
 603                   }
 604
 605                 }
 606               }
 607             }//end loop layer 6
 608           }//end loop layer 5
 609         }//end loop layer 4
 610       }//end loop layer 3
 611     }//end loop layer 2
 612   }//end loop layer 1
 613
 614
 615
 616
 617   if(fListOfTracks->GetEntries()==0) return 0;
 618
 619   Int_t lowchi2 = FindTrackLowChiSquare();
 620   AliITStrackV2* otrack =(AliITStrackV2*)fListOfTracks->At(lowchi2);
 621   AliITStrackSA* trsa = (AliITStrackSA*)fListOfSATracks->At(lowchi2);
 622
 623   if(otrack==0) return 0;
 624
 625   CookLabel(otrack,0.); //MI change - to see fake ratio
 626   Int_t label=FindLabel(otrack);
 627   otrack->SetLabel(label);
 628   Double_t low=0.;
 629   Double_t up=0.51;
 630   otrack->CookdEdx(low,up);
 631
 632   //remove clusters of found track
 633   for(Int_t nlay=0;nlay<AliITSgeomTGeo::GetNLayers();nlay++){
 634     for(Int_t cln=0;cln<trsa->GetNumberOfMarked(nlay);cln++){
 635       Int_t index = trsa->GetClusterMark(nlay,cln);
 636       RemoveClusterCoord(nlay,index);
 637     }
 638   }
 639
 640   return otrack;
 641
 642 }
 643
 644 //_______________________________________________________
 645 void AliITStrackerSA::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) const
 646 {
 647   //
 648   // Add new track to the ESD
 649   //
 650   AliESDtrack outtrack;
 651   outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
 652   if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
 653   for(Int_t i=0;i<12;i++) {
 654     outtrack.SetITSModuleIndex(i,t->GetModuleIndex(i));
 655   }
 656   Double_t sdedx[4]={0.,0.,0.,0.};
 657   for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
 658   outtrack.SetITSdEdxSamples(sdedx);
 659
 660
 661   if(AliITSReconstructor::GetRecoParam()->GetSAUsedEdxInfo()){
 662     Double_t mom=t->P();
 663     Double_t ppid[AliPID::kSPECIES];
 664     for(Int_t isp=0;isp<AliPID::kSPECIES;isp++) ppid[isp]=0.;
 665     ppid[AliPID::kPion]=1.;
 666     if(mom<0.7){
 667       Double_t truncmean=t->GetdEdx();
 668       Int_t ide=fITSPid->GetParticleIdFromdEdxVsP(mom,truncmean,kTRUE);
 669       if(ide==AliPID::kProton){
 670         ppid[AliPID::kProton]=1.;
 671         ppid[AliPID::kPion]=0.;
 672       }
 673       else if(ide==AliPID::kKaon){
 674         ppid[AliPID::kKaon]=1.;
 675         ppid[AliPID::kPion]=0.;
 676       }
 677     }
 678     outtrack.SetITSpid(ppid);
 679     outtrack.SetESDpid(ppid);
 680   }
 681   event->AddTrack(&outtrack);
 682
 683   return;
 684 }
 685
 686
 687 //_______________________________________________________
 688 Int_t AliITStrackerSA::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackSA* trs,Double_t /*zvertex*/,Int_t pflag){
 689   //function used to to find the clusters associated to the track
 690
 691   if(ForceSkippingOfLayer(layer)) return 0;
 692
 693   Int_t nc=0;
 694   AliITSlayer &lay = fgLayers[layer];
 695   Double_t r=lay.GetR();
 696   if(pflag==1){
 697     Double_t cx1,cx2,cy1,cy2;
 698     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
 699     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
 700        return 0;
 701     Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
 702     Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
 703     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
 704   }
 705
 706
 707   Double_t phiExpect=fPhiEstimate;
 708   Double_t lamExpect=fLambdac;
 709
 710   Int_t ncl = fCluCoord[layer]->GetEntriesFast();
 711   for (Int_t index=0; index<ncl; index++) {
 712
 713     AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
 714
 715     Double_t lambda = arr->GetLambda();
 716     if (TMath::Abs(lambda-lamExpect)>lambdawindow) continue;
 717
 718     Double_t phi = arr->GetPhi();
 719     Double_t deltaPhi = phi-phiExpect;
 720     if(deltaPhi>TMath::Pi()) deltaPhi-=2*TMath::Pi();
 721     else if(deltaPhi<-TMath::Pi()) deltaPhi+=2*TMath::Pi();
 722     if (TMath::Abs(deltaPhi)>phiwindow) continue;
 723
 724     if(trs->GetNumberOfClustersSA()==trs->GetMaxNumberOfClusters()) return 0;
 725     if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
 726     Int_t orind = arr->GetOrInd();
 727     trs->AddClusterSA(layer,orind);
 728     trs->AddClusterMark(layer,index);
 729     nc++;
 730     fLambdac=lambda;
 731     fPhiEstimate=phi;
 732
 733     fPointc[0]=arr->GetX();
 734     fPointc[1]=arr->GetY();
 735
 736   }
 737   return nc;
 738 }
 739
 740 //________________________________________________________________
 741 Bool_t AliITStrackerSA::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
 742   // Sets the first point (seed) for tracking
 743
 744   AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
 745   fPhic = arr->GetPhi();
 746   fLambdac = arr->GetLambda();
 747   fPhiEstimate = fPhic;
 748   fPoint1[0]=primaryVertex[0];
 749   fPoint1[1]=primaryVertex[1];
 750   fPoint2[0]=arr->GetX();
 751   fPoint2[1]=arr->GetY();
 752   return kTRUE;
 753 }
 754
 755 //________________________________________________________________
 756 void AliITStrackerSA::UpdatePoints(){
 757   //update of points for the estimation of the curvature
 758
 759   fPoint2[0]=fPoint3[0];
 760   fPoint2[1]=fPoint3[1];
 761   fPoint3[0]=fPointc[0];
 762   fPoint3[1]=fPointc[1];
 763
 764
 765 }
 766
 767 //___________________________________________________________________
 768 Int_t AliITStrackerSA::FindEquation(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3,Double_t& a, Double_t& b, Double_t& c){
 769
 770    //given (x,y) of three recpoints (in global coordinates)
 771    //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
 772
 773    Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 774    if(den==0) return 0;
 775    a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 776    b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 777    c = -x1*x1-y1*y1-a*x1-b*y1;
 778    return 1;
 779  }
 780 //__________________________________________________________________________
 781  Int_t AliITStrackerSA::FindIntersection(Double_t a1, Double_t b1, Double_t c1, Double_t c2,Double_t& x1,Double_t& y1, Double_t& x2, Double_t& y2){
 782
 783  //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
 784  //c2 is -rlayer*rlayer
 785
 786   if(a1==0) return 0;
 787  Double_t m = c2-c1;
 788  Double_t aA = (b1*b1)/(a1*a1)+1;
 789  Double_t bB = (-2*m*b1/(a1*a1));
 790  Double_t cC = c2+(m*m)/(a1*a1);
 791  Double_t dD = bB*bB-4*aA*cC;
 792  if(dD<0) return 0;
 793
 794  y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
 795  y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
 796  x1 = (c2-c1-b1*y1)/a1;
 797  x2 = (c2-c1-b1*y2)/a1;
 798
 799  return 1;
 800 }
 801 //____________________________________________________________________
 802 Double_t AliITStrackerSA::Curvature(Double_t x1,Double_t y1,Double_t
 803 x2,Double_t y2,Double_t x3,Double_t y3){
 804
 805   //calculates the curvature of track
 806   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 807   if(den==0) return 0;
 808   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 809   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 810   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
 811   Double_t xc=-a/2.;
 812
 813   if((a*a+b*b-4*c)<0) return 0;
 814   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
 815   if(rad==0) return 0;
 816
 817   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
 818   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
 819   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
 820   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
 821
 822   return 1/rad;
 823
 824 }
 825
 826
 827 //____________________________________________________________________
 828 Double_t AliITStrackerSA::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
 829
 830   //Returns the point closest to pp
 831
 832   Double_t diff1 = p1-pp;
 833   Double_t diff2 = p2-pp;
 834
 835   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
 836   else fPhiEstimate=p2;
 837   return fPhiEstimate;
 838
 839 }
 840
 841
 842 //_________________________________________________________________
 843 Int_t AliITStrackerSA::FindTrackLowChiSquare() const {
 844   // returns track with lowest chi square
 845   Int_t dim=fListOfTracks->GetEntries();
 846   if(dim<=1) return 0;
 847   AliITStrackV2* trk = (AliITStrackV2*)fListOfTracks->At(0);
 848   Double_t minChi2=trk->GetChi2();
 849   Int_t index=0;
 850   for(Int_t i=1;i<dim;i++){
 851     trk = (AliITStrackV2*)fListOfTracks->At(i);
 852     Double_t chi2=trk->GetChi2();
 853     if(chi2<minChi2){
 854       minChi2=chi2;
 855       index=i;
 856     }
 857   }
 858   return index;
 859 }
 860
 861 //__________________________________________________________
 862 Int_t AliITStrackerSA::FindLabel(AliITStrackV2* track){
 863   // compute the track label starting from cluster labels
 864
 865   Int_t labl[AliITSgeomTGeo::kNLayers][3];
 866   Int_t cnts[AliITSgeomTGeo::kNLayers][3];
 867   for(Int_t j=0;j<AliITSgeomTGeo::GetNLayers();j++){
 868     for(Int_t k=0;k<3;k++){
 869       labl[j][k]=-2;
 870       cnts[j][k]=1;
 871     }
 872   }
 873   Int_t iNotLabel=0;
 874   for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
 875     Int_t indexc = track->GetClusterIndex(i);
 876     AliITSRecPoint* cl = (AliITSRecPoint*)GetCluster(indexc);
 877     Int_t iLayer=cl->GetLayer();
 878     for(Int_t k=0;k<3;k++){
 879       labl[iLayer][k]=cl->GetLabel(k);
 880       if(labl[iLayer][k]<0) iNotLabel++;
 881     }
 882   }
 883   if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
 884
 885   for(Int_t j1=0;j1<AliITSgeomTGeo::kNLayers; j1++) {
 886     for(Int_t j2=0; j2<j1;  j2++){
 887       for(Int_t k1=0; k1<3; k1++){
 888         for(Int_t k2=0; k2<3; k2++){
 889           if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
 890             cnts[j2][k2]++;
 891             cnts[j1][k1]=0;
 892           }
 893         }
 894       }
 895     }
 896   }
 897
 898
 899   Int_t cntMax=0;
 900   Int_t label=-1;
 901   for(Int_t j=0;j<AliITSgeomTGeo::kNLayers;j++){
 902     for(Int_t k=0;k<3;k++){
 903       if(cnts[j][k]>cntMax && labl[j][k]>=0){
 904         cntMax=cnts[j][k];
 905         label=labl[j][k];
 906       }
 907     }
 908   }
 909
 910   Int_t lflag=0;
 911   for(Int_t i=0;i<AliITSgeomTGeo::kNLayers;i++)
 912     if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
 913
 914   if(lflag<track->GetNumberOfClusters()) label = -label;
 915   return label;
 916 }
 917 //_____________________________________________________________________________
 918 void AliITStrackerSA::SetCalculatedWindowSizes(Int_t n, Double_t phimin, Double_t phimax, Double_t lambdamin, Double_t lambdamax){
 919   // Set sizes of the phi and lambda windows used for track finding
 920   fNloop = n;
 921   if(fPhiWin) delete [] fPhiWin;
 922   if(fLambdaWin) delete [] fLambdaWin;
 923   fPhiWin = new Double_t[fNloop];
 924   fLambdaWin = new Double_t[fNloop];
 925   Double_t stepPhi=(phimax-phimin)/(Double_t)(fNloop-1);
 926   Double_t stepLambda=(lambdamax-lambdamin)/(Double_t)(fNloop-1);
 927   for(Int_t k=0;k<fNloop;k++){
 928     Double_t phi=phimin+k*stepPhi;
 929     Double_t lam=lambdamin+k*stepLambda;
 930     fPhiWin[k]=phi;
 931     fLambdaWin[k]=lam;
 932   }
 933 }
 934 //_____________________________________________________________________________
 935 void AliITStrackerSA::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
 936   // Set sizes of the phi and lambda windows used for track finding
 937   fNloop = n;
 938   if(phi){ // user defined values
 939     fPhiWin = new Double_t[fNloop];
 940     fLambdaWin = new Double_t[fNloop];
 941     for(Int_t k=0;k<fNloop;k++){
 942       fPhiWin[k]=phi[k];
 943       fLambdaWin[k]=lam[k];
 944     }
 945   }
 946   else {  // default values
 947
 948     Double_t phid[33]   = {0.002,0.003,0.004,0.0045,0.0047,
 949                            0.005,0.0053,0.0055,
 950                            0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
 951                            0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
 952                            0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
 953     Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
 954                             0.005,0.005,0.006,
 955                             0.006,0.006,0.006,0.007,0.007,0.007,0.007,
 956                             0.007,0.007,0.007,0.007,0.007,0.007,0.007,
 957                             0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
 958
 959     if(fNloop!=33){
 960       fNloop = 33;
 961     }
 962
 963
 964     fPhiWin = new Double_t[fNloop];
 965     fLambdaWin = new Double_t[fNloop];
 966
 967     Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes(); // possibility to enlarge windows for cosmics reco with large misalignments (A.Dainese)
 968
 969     for(Int_t k=0;k<fNloop;k++){
 970       fPhiWin[k]=phid[k]*factor;
 971       fLambdaWin[k]=lambdad[k]*factor;
 972     }
 973
 974   }
 975
 976 }
 977 //_______________________________________________________________________
 978 void AliITStrackerSA::GetCoorAngles(AliITSRecPoint* cl,Double_t &phi,Double_t &lambda, Double_t &x, Double_t &y,Double_t &z, const Double_t* vertex){
 979   //Returns values of phi (azimuthal) and lambda angles for a given cluster
 980 /*
 981   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
 982   Int_t lay,lad,det; fGeom->GetModuleId(module,lay,lad,det);
 983   Double_t tx,ty,tz;  fGeom->GetTrans(lay,lad,det,tx,ty,tz);
 984
 985   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
 986   Double_t phi1=TMath::Pi()/2+alpha;
 987   if (lay==1) phi1+=TMath::Pi();
 988
 989   Double_t cp=TMath::Cos(phi1), sp=TMath::Sin(phi1);
 990   Double_t r=tx*cp+ty*sp;
 991
 992   xyz= r*cp - cl->GetY()*sp;
 993   y= r*sp + cl->GetY()*cp;
 994   z=cl->GetZ();
 995 */
 996   Float_t xyz[3];
 997   cl->GetGlobalXYZ(xyz);
 998   x=xyz[0];
 999   y=xyz[1];
1000   z=xyz[2];
1001
1002   phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1003   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1004 }
1005
1006 //________________________________________________________________________
1007 void AliITStrackerSA::GetCoorErrors(AliITSRecPoint* cl,Double_t &sx,Double_t &sy, Double_t &sz){
1008
1009   //returns sigmax, y, z of cluster in global coordinates
1010 /*
1011   Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1012   Int_t lay,lad,det;
1013   AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1014
1015   Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1016   Double_t phi=TMath::Pi()/2+alpha;
1017   if (lay==1) phi+=TMath::Pi();
1018
1019   Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1020 */
1021   Float_t covm[6];
1022   cl->GetGlobalCov(covm);
1023   sx=TMath::Sqrt(covm[0]);
1024   sy=TMath::Sqrt(covm[3]);
1025   sz=TMath::Sqrt(covm[5]);
1026 /*
1027   sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1028   sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1029   sz = TMath::Sqrt(cl->GetSigmaZ2());
1030 */
1031 }
1032