ITS/UPGRADE/AliITStrackerUpgrade.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 //  ITS Upgrade Stand alone tracker class         //
  20 //  Authors: A.Mastroserio C.Terrevoli            //
  21 //  e-mail:annalisa.mastroserio@cern.ch           //
  22 //         cristina.terrevoli@ba.infn.it          //
  23 ////////////////////////////////////////////////////
  24
  25 #include <stdlib.h>
  26
  27 #include <TArrayI.h>
  28 #include <TBranch.h>
  29 #include <TObjArray.h>
  30 #include <TTree.h>
  31 #include "AliRunLoader.h"
  32 #include "AliITSLoader.h"
  33 #include "AliESDEvent.h"
  34 #include "AliESDVertex.h"
  35 #include "AliESDtrack.h"
  36 #include "AliITSVertexer.h"
  37 #include "AliITSclusterTable.h"
  38 #include "AliITSRecPointU.h"
  39 #include "AliITStrackU.h"
  40 #include "AliITStrackerMI.h"
  41 #include "AliITSlayerUpgrade.h"
  42 #include "AliITSsegmentationUpgrade.h"
  43 #include "AliITStrackerUpgrade.h"
  44 #include "AliITSReconstructor.h"
  45 #include "AliLog.h"
  46 #include "AliRun.h"
  47
  48 ClassImp(AliITStrackerUpgrade)
  49
  50 //____________________________________________________________________________
  51   AliITStrackerUpgrade::AliITStrackerUpgrade():AliITStrackerMI(),
  52                                                fNLayers(), fPhiEstimate(0), fITSStandAlone(0), fLambdac(0),
  53                                                fPhic(0), fCoef1(0), fCoef2(0), fCoef3(0), fNloop(0),
  54                                                fPhiWin(0),  fLambdaWin(0), fVert(0), fVertexer(0),
  55                                                fListOfTracks(0), fListOfUTracks(0), fITSclusters(0),
  56                                                fInwardFlag(0), fOuterStartLayer(0),  fInnerStartLayer(0),
  57                                                fMinNPoints(0),  fMinQ(0), fTrackToFollow(), fLayers(0), fSegmentation(0x0),
  58                                                fCluLayer(0), fCluCoord(0)
  59 {
  60   // Default constructor
  61   Init();
  62
  63 }
  64 //____________________________________________________________________________
  65 AliITStrackerUpgrade::AliITStrackerUpgrade(Int_t nLay):AliITStrackerMI(),
  66                                                        fNLayers(nLay), fPhiEstimate(0), fITSStandAlone(0), fLambdac(0),
  67                                                        fPhic(0), fCoef1(0), fCoef2(0), fCoef3(0), fNloop(0),
  68                                                        fPhiWin(0),  fLambdaWin(0), fVert(0), fVertexer(0),
  69                                                        fListOfTracks(0), fListOfUTracks(0), fITSclusters(0),
  70                                                        fInwardFlag(0), fOuterStartLayer(0),  fInnerStartLayer(nLay),
  71                                                        fMinNPoints(0),  fMinQ(0),  fTrackToFollow(), fLayers(0), fSegmentation(0x0),
  72                                                        fCluLayer(0), fCluCoord(0)
  73 {
  74   //
  75   // constructor
  76   //
  77   Init();
  78
  79 }
  80
  81 //________________________________________________________________________
  82 AliITStrackerUpgrade::AliITStrackerUpgrade(const AliITStrackerUpgrade& tracker):AliITStrackerMI(),
  83                                                                                 fNLayers(tracker.fNLayers), fPhiEstimate(tracker.fPhiEstimate),
  84                                                                                 fITSStandAlone(tracker.fITSStandAlone), fLambdac(tracker.fLambdac),
  85                                                                                 fPhic(tracker.fPhic), fCoef1(tracker.fCoef1), fCoef2(tracker.fCoef2),
  86                                                                                 fCoef3(tracker.fCoef3), fNloop(tracker.fNloop), fPhiWin(tracker.fPhiWin),
  87                                                                                 fLambdaWin(tracker.fLambdaWin), fVert(tracker.fVert), fVertexer(tracker.fVertexer),
  88                                                                                 fListOfTracks(tracker.fListOfTracks), fListOfUTracks(tracker.fListOfUTracks),
  89                                                                                 fITSclusters(tracker.fITSclusters), fInwardFlag(tracker.fInwardFlag),
  90                                                                                 fOuterStartLayer(tracker.fOuterStartLayer), fInnerStartLayer(tracker.fInnerStartLayer),
  91                                                                                 fMinNPoints(tracker.fMinNPoints), fMinQ(tracker.fMinQ), fTrackToFollow(tracker.fTrackToFollow), fLayers(tracker.fLayers),
  92                                                                                 fSegmentation(tracker.fSegmentation), fCluLayer(tracker.fCluLayer), fCluCoord(tracker.fCluCoord)
  93 {
  94   //
  95   // Copy constructor
  96   for(Int_t i=0;i<2;i++){
  97     fPoint1[i]=tracker.fPoint1[i];
  98     fPoint2[i]=tracker.fPoint2[i];
  99     fPoint3[i]=tracker.fPoint3[i];
 100     fPointc[i]=tracker.fPointc[i];
 101   }
 102   if(tracker.fVertexer && tracker.fVert){
 103     fVert = new AliESDVertex(*tracker.fVert);
 104   }
 105   else {
 106     fVert = tracker.fVert;
 107   }
 108   for(Int_t i=0;i<fNLayers;i++){
 109     fCluLayer[i] = tracker.fCluLayer[i];
 110     fCluCoord[i] = tracker.fCluCoord[i];
 111   }
 112 }
 113 //______________________________________________________________________
 114 AliITStrackerUpgrade& AliITStrackerUpgrade::operator=(const AliITStrackerUpgrade& source){
 115   // Assignment operator.
 116   this->~AliITStrackerUpgrade();
 117   new(this) AliITStrackerUpgrade(source);
 118   return *this;
 119
 120 }
 121
 122 //____________________________________________________________________________
 123 AliITStrackerUpgrade::~AliITStrackerUpgrade(){
 124   // destructor
 125   // if fVertexer is not null, the AliESDVertex obj. is owned by this class
 126   // and is deleted here
 127   if(fVertexer){
 128     if(fVert)delete fVert;
 129   }
 130   fVert = 0;
 131   fVertexer = 0;
 132
 133   if(fPhiWin)delete []fPhiWin;
 134   if(fLambdaWin)delete []fLambdaWin;
 135   fListOfTracks->Delete();
 136   delete fListOfTracks;
 137   fListOfUTracks->Delete();
 138   delete fListOfUTracks;
 139   if(fCluLayer){
 140     for(Int_t i=0;i<fNLayers;i++){
 141       if(fCluLayer[i]){
 142         fCluLayer[i]->Delete();
 143         delete fCluLayer[i];
 144       }
 145     }
 146     delete [] fCluLayer;
 147   }
 148   if(fCluCoord){
 149     for(Int_t i=0;i<fNLayers;i++){
 150       if(fCluCoord[i]){
 151         fCluCoord[i]->Delete();
 152         delete fCluCoord[i];
 153       }
 154     }
 155     delete [] fCluCoord;
 156   }
 157
 158   if(fSegmentation) delete fSegmentation;
 159 }
 160
 161 //____________________________________________________________________________
 162 Int_t AliITStrackerUpgrade::Clusters2Tracks(AliESDEvent *event){
 163   // This method is used to find and fit the tracks. By default the corresponding
 164   // method in the parent class is invoked. In this way a combined tracking
 165   // TPC+ITS is performed. If the flag fITSStandAlone is true, the tracking
 166   // is done in the ITS only. In the standard reconstruction chain this option
 167   // can be set via AliReconstruction::SetOption("ITS","onlyITS")
 168   Int_t rc=0;
 169   if(!fITSStandAlone){
 170     rc=AliITStrackerMI::Clusters2Tracks(event);
 171   }
 172   else {
 173     AliDebug(1,"Stand Alone flag set: doing tracking in ITS alone\n");
 174   }
 175   if(!rc){
 176     rc=FindTracks(event,kFALSE);
 177     if(AliITSReconstructor::GetRecoParam()->GetSAUseAllClusters()==kTRUE) {
 178       rc=FindTracks(event,kTRUE);
 179     }
 180   }
 181   return rc;
 182 }
 183 //_________________________________________________-
 184 void AliITStrackerUpgrade::Init(){
 185   //  Reset all data members
 186   const Int_t nL = fNLayers;
 187   fPhiEstimate=0;
 188   for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 189   fLambdac=0;
 190   fPhic=0;
 191   fCoef1=0;
 192   fCoef2=0;
 193   fCoef3=0;
 194   fPointc[0]=0;
 195   fPointc[1]=0;
 196   fVert = 0;
 197   fVertexer = 0;
 198   Int_t nLoops=AliITSReconstructor::GetRecoParam()->GetNLoopsSA();
 199   if(nLoops==33){
 200     SetFixedWindowSizes();
 201   }else{
 202     Double_t phimin=AliITSReconstructor::GetRecoParam()->GetMinPhiSA();
 203     Double_t phimax=AliITSReconstructor::GetRecoParam()->GetMaxPhiSA();
 204     Double_t lambmin=AliITSReconstructor::GetRecoParam()->GetMinLambdaSA();
 205     Double_t lambmax=AliITSReconstructor::GetRecoParam()->GetMaxLambdaSA();
 206     SetCalculatedWindowSizes(nLoops,phimin,phimax,lambmin,lambmax);
 207   }
 208   fMinQ=AliITSReconstructor::GetRecoParam()->GetSAMinClusterCharge();
 209   fITSclusters = 0;
 210   SetOuterStartLayer(1);
 211   SetSAFlag(kFALSE);
 212   fListOfTracks=new TClonesArray("AliITStrackU",100);
 213   fListOfUTracks=new TClonesArray("AliITStrackU",100);
 214   fLayers=new AliITSlayerUpgrade*[nL];//to be fixed
 215   fCluLayer = 0;
 216   fCluCoord = 0;
 217   fMinNPoints = 3;
 218   for(Int_t layer=0; layer<fNLayers; layer++){
 219     Double_t p=0.;
 220     Double_t zC= 0.;
 221     fLayers[layer] = new AliITSlayerUpgrade(p,zC);
 222   }
 223
 224   fSegmentation = new AliITSsegmentationUpgrade();
 225
 226 }
 227 //_______________________________________________________________________
 228 Int_t AliITStrackerUpgrade::LoadClusters(TTree *clusTree){
 229   //
 230   // Load clusters for tracking
 231   //
 232
 233   TClonesArray statITSCluster("AliITSRecPointU");
 234   TClonesArray *ITSCluster = &statITSCluster;
 235
 236   TBranch* itsClusterBranch=clusTree->GetBranch("ITSRecPoints");
 237   if (!itsClusterBranch){
 238     AliError("can't get the branch with the ITS clusters ! \n");
 239     return 1;
 240   }
 241   itsClusterBranch->SetAddress(&ITSCluster);
 242   clusTree->GetEvent(0);
 243   Int_t nCluster = ITSCluster->GetEntriesFast();
 244   for(Int_t i=0; i<nCluster; i++){
 245     AliITSRecPointU *recp = (AliITSRecPointU*)ITSCluster->UncheckedAt(i);
 246     fLayers[recp->GetLayer()]->InsertCluster(new AliITSRecPointU(*recp));
 247   }//loop clusters
 248
 249   SetClusterTree(clusTree);
 250   return 0;
 251 }
 252
 253 //_______________________________________________________
 254 void AliITStrackerUpgrade::ResetForFinding(){
 255   //  Reset data members used in all loops during track finding
 256   fPhiEstimate=0;
 257   for(Int_t i=0;i<3;i++){fPoint1[i]=0;fPoint2[i]=0;fPoint3[i]=0;}
 258   fLambdac=0;
 259   fPhic=0;
 260   fCoef1=0;
 261   fCoef2=0;
 262   fCoef3=0;
 263   fPointc[0]=0;
 264   fPointc[1]=0;
 265   fListOfTracks->Clear();
 266   fListOfUTracks->Clear();
 267 }
 268
 269
 270
 271 //______________________________________________________________________
 272 Int_t AliITStrackerUpgrade::FindTracks(AliESDEvent* event,Bool_t useAllClusters){
 273
 274   // Track finder using the ESD object
 275   AliDebug(2,Form(" field is %f",event->GetMagneticField()));
 276   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
 277   if(!fITSclusters){
 278     Fatal("FindTracks","ITS cluster tree is not accessed - Abort!!!\n Please use method SetClusterTree to pass the pointer to the tree\n");
 279     return -1;
 280   }
 281
 282   //Reads event and mark clusters of traks already found, with flag kITSin
 283   Int_t nentr=event->GetNumberOfTracks();
 284   if(!useAllClusters) {
 285     while (nentr--) {
 286       AliESDtrack *track=event->GetTrack(nentr);
 287       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
 288         Int_t idx[12];
 289         Int_t ncl = track->GetITSclusters(idx);
 290         for(Int_t k=0;k<ncl;k++){
 291           AliITSRecPointU* cll = (AliITSRecPointU*)GetCluster(idx[k]);
 292           cll->SetBit(kSAflag);
 293         }
 294       }
 295     }
 296   }else{
 297     while (nentr--) {
 298       AliESDtrack *track=event->GetTrack(nentr);
 299       if ((track->GetStatus()&AliESDtrack::kITSin) == AliESDtrack::kITSin){
 300         Int_t idx[12];
 301         Int_t ncl = track->GetITSclusters(idx);
 302         for(Int_t k=0;k<ncl;k++){
 303           AliITSRecPointU* cll = (AliITSRecPointU*)GetCluster(idx[k]);
 304           cll->ResetBit(kSAflag);
 305         }
 306       }
 307     }
 308   }
 309   //Get primary vertex
 310   Double_t primaryVertex[3];
 311   event->GetVertex()->GetXYZ(primaryVertex);
 312   Int_t nclusters[fNLayers];
 313   for(Int_t i=0; i<fNLayers; i++){
 314     nclusters[i]=0;
 315   }
 316
 317   Int_t dmar[fNLayers];
 318   for(Int_t i=0; i<fNLayers; i++){
 319     dmar[i]=0;
 320   }
 321
 322   if (fCluLayer == 0) {
 323     fCluLayer = new TClonesArray*[fNLayers];
 324     fCluCoord = new TClonesArray*[fNLayers];
 325     for(Int_t i=0;i<fNLayers;i++) {
 326       fCluLayer[i]=0;
 327       fCluCoord[i]=0;
 328     }
 329   }
 330   for(Int_t i=0;i<fNLayers;i++){
 331     if (!ForceSkippingOfLayer(i)) {
 332       for(Int_t cli=0;cli<fLayers[i]->GetNumberOfClusters();cli++){
 333         AliITSRecPointU* cls = (AliITSRecPointU*)fLayers[i]->GetCluster(cli);
 334         if(cls->TestBit(kSAflag)==kTRUE) continue; //clusters used by TPC prol.
 335         if(cls->GetQ()==0) continue; //fake clusters dead zones
 336         if(i>1 && cls->GetQ()<=fMinQ) continue; // cut on SDD and SSD cluster charge
 337         nclusters[i]++;
 338       }
 339     }
 340     dmar[i]=0;
 341     if(!fCluLayer[i]){
 342       fCluLayer[i] = new TClonesArray("AliITSRecPointU",nclusters[i]);
 343     }else{
 344       fCluLayer[i]->Delete();
 345       fCluLayer[i]->Expand(nclusters[i]);
 346     }
 347     if(!fCluCoord[i]){
 348       fCluCoord[i] = new TClonesArray("AliITSclusterTable",nclusters[i]);
 349     }else{
 350       fCluCoord[i]->Delete();
 351       fCluCoord[i]->Expand(nclusters[i]);
 352     }
 353   }
 354   for(Int_t ilay=0;ilay<fNLayers;ilay++){
 355     TClonesArray &clulay = *fCluLayer[ilay];
 356     TClonesArray &clucoo = *fCluCoord[ilay];
 357     if (!ForceSkippingOfLayer(ilay)){
 358       AliDebug(2,Form("number of clusters in layer %i : %i",ilay,fLayers[ilay]->GetNumberOfClusters()));
 359       for(Int_t cli=0;cli<fLayers[ilay]->GetNumberOfClusters();cli++){
 360         AliITSRecPointU* cls = (AliITSRecPointU*)fLayers[ilay]->GetCluster(cli);
 361         if(cls->TestBit(kSAflag)==kTRUE) continue;
 362         if(cls->GetQ()==0) continue;
 363         Double_t phi=0;Double_t lambda=0;
 364         Double_t x=0;Double_t y=0;Double_t z=0;
 365         Float_t sx=0;Float_t sy=0;Float_t sz=0;
 366         GetCoorAngles(cls,phi,lambda,x,y,z,primaryVertex);
 367         GetCoorErrors(cls,sx,sy,sz);
 368         new (clulay[dmar[ilay]]) AliITSRecPoint(*cls);
 369         new (clucoo[dmar[ilay]]) AliITSclusterTable(x,y,z,sx,sy,sz,phi,lambda,cli);
 370         dmar[ilay]++;
 371       }
 372     } else AliDebug(2,Form("Force skipping layer %i",ilay));
 373   }
 374
 375
 376   // track counter
 377   Int_t ntrack=0;
 378
 379   static Int_t nClusLay[fgMaxNLayer];//counter for clusters on each layer
 380   Int_t startLayForSeed=0;
 381   Int_t lastLayForSeed=fOuterStartLayer;
 382   Int_t nSeedSteps=lastLayForSeed-startLayForSeed;
 383   Int_t seedStep=1;
 384   if(fInwardFlag){
 385     startLayForSeed=fNLayers-1;
 386     lastLayForSeed=fInnerStartLayer;
 387     nSeedSteps=startLayForSeed-lastLayForSeed;
 388     seedStep=-1;
 389   }
 390   // loop on minimum number of points
 391   for(Int_t iMinNPoints=fNLayers; iMinNPoints>=fMinNPoints; iMinNPoints--) {
 392     // loop on starting layer for track finding
 393     for(Int_t iSeedLay=0; iSeedLay<=nSeedSteps; iSeedLay++) {
 394       Int_t theLay=startLayForSeed+iSeedLay*seedStep;
 395       if(ForceSkippingOfLayer(theLay)) continue;
 396       Int_t minNPoints=iMinNPoints-theLay;
 397       if(fInwardFlag) minNPoints=iMinNPoints-(fNLayers-1-theLay);
 398       for(Int_t i=theLay+1;i<fNLayers;i++)
 399         if(ForceSkippingOfLayer(i))
 400           minNPoints--;
 401       if(minNPoints<fMinNPoints) continue;
 402
 403       // loop on phi and lambda window size
 404       for(Int_t nloop=0;nloop<fNloop;nloop++){
 405         Int_t nclTheLay=fCluLayer[theLay]->GetEntries();
 406         while(nclTheLay--){
 407           ResetForFinding();
 408           Bool_t useRP=SetFirstPoint(theLay,nclTheLay,primaryVertex);
 409           if(!useRP) continue;
 410           AliITStrackU trs(fNLayers);
 411
 412           Int_t pflag=0;
 413           Int_t kk;
 414           for(kk=0;kk<fNLayers;kk++) nClusLay[kk] = 0;
 415           kk=0;
 416           nClusLay[kk] = SearchClusters(theLay,fPhiWin[nloop],fLambdaWin[nloop],
 417                                         &trs,primaryVertex[2],pflag);
 418           Int_t nextLay=theLay+seedStep;
 419           Bool_t goon=kTRUE;
 420           while(goon){
 421
 422             kk++;
 423             nClusLay[kk] = SearchClusters(nextLay,fPhiWin[nloop],fLambdaWin[nloop],
 424                                           &trs,primaryVertex[2],pflag);
 425             if(nClusLay[kk]!=0){
 426
 427               pflag=1;
 428               if(kk==1) {
 429                 fPoint3[0]=fPointc[0];
 430                 fPoint3[1]=fPointc[1];
 431               } else {
 432                 UpdatePoints();
 433               }
 434             }
 435             nextLay+=seedStep;
 436             if(nextLay<0 || nextLay==fNLayers) goon=kFALSE;
 437           }
 438
 439
 440           Int_t layOK=0;
 441           if(!fInwardFlag){
 442             for(Int_t nnp=0;nnp<fNLayers-theLay;nnp++){
 443               if(nClusLay[nnp]!=0) layOK+=1;
 444             }
 445           }else{
 446             for(Int_t nnp=theLay; nnp>=0; nnp--){
 447               if(nClusLay[nnp]!=0) layOK+=1;
 448             }
 449           }
 450           if(layOK>=minNPoints){
 451             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]));
 452             AliITStrackV2* tr2 = 0;
 453             Bool_t onePoint = kFALSE;
 454             tr2 = FitTrack(&trs,primaryVertex,onePoint);
 455             if(!tr2){
 456               continue;
 457             }
 458             AliDebug(2,Form("---NPOINTS fit: %d\n",tr2->GetNumberOfClusters()));
 459
 460             StoreTrack(tr2,event,useAllClusters);
 461             ntrack++;
 462
 463           }
 464
 465         }//end loop on clusters of theLay
 466       } //end loop on window sizes
 467     } //end loop on theLay
 468   }//end loop on min points
 469   if(!useAllClusters) AliInfo(Form("Number of found tracks: %d",event->GetNumberOfTracks()));
 470
 471   ResetForFinding();
 472   return 0;
 473
 474 }
 475
 476 //________________________________________________________________________
 477
 478 AliITStrackV2* AliITStrackerUpgrade::FitTrack(AliITStrackU* tr,Double_t *primaryVertex,Bool_t onePoint) {
 479
 480   const Int_t kMaxClu=AliITStrackU::kMaxNumberOfClusters;
 481
 482
 483   static Int_t clind[fgMaxNLayer][kMaxClu];
 484   static Int_t clmark[fgMaxNLayer][kMaxClu];
 485   static Int_t end[fgMaxNLayer];
 486   static AliITSRecPointU *listlayer[fgMaxNLayer][kMaxClu];
 487   for(Int_t k=0;k<fgMaxNLayer; k++)end[k]=0;
 488
 489   for(Int_t i=0;i<fNLayers;i++) {
 490     for(Int_t j=0;j<kMaxClu; j++){
 491       clind[i][j]=0;
 492       clmark[i][j]=0;
 493       listlayer[i][j]=0;
 494     }
 495   }
 496   Int_t inx[fNLayers];
 497   for (Int_t k=0; k<fNLayers; k++) inx[k]=-1;
 498   Int_t nclusters = tr->GetNumberOfClustersU();
 499   for(Int_t ncl=0;ncl<nclusters;ncl++){
 500     Int_t index = tr->GetClusterIndexU(ncl);
 501     Int_t lay = (index & 0xf0000000) >> 28;
 502     //Int_t cli = index&0x0fffffff;
 503     //AliITSRecPoint* cl = (AliITSRecPoint*)fLayers[lay]->GetCluster(cli);
 504     AliITSRecPointU* cl = (AliITSRecPointU*)GetCluster(index);
 505     if(cl->TestBit(kSAflag)==kTRUE) cl->ResetBit(kSAflag);
 506     inx[lay]=index;
 507
 508     Int_t nInLay=end[lay];
 509     listlayer[lay][nInLay]=cl;
 510     clind[lay][nInLay]=index;
 511     end[lay]++;
 512   }
 513
 514   for(Int_t nlay=0;nlay<fNLayers;nlay++){
 515     for(Int_t ncl=0;ncl<tr->GetNumberOfMarked(nlay);ncl++){
 516       Int_t mark = tr->GetClusterMark(nlay,ncl);
 517       clmark[nlay][ncl]=mark;
 518     }
 519   }
 520
 521
 522   Int_t firstLay=-1,secondLay=-1;
 523   for(Int_t i=0;i<fNLayers;i++) {
 524     if(end[i]==0) {
 525       end[i]=1;
 526     }else{
 527       if(firstLay==-1) {
 528         firstLay=i;
 529       } else if(secondLay==-1) {
 530         secondLay=i;
 531       }
 532     }
 533   }
 534
 535   //---adding initialization for non existent layers
 536   for(Int_t i=fNLayers;i<fgMaxNLayer;i++) {
 537     if(end[i]==0) {
 538       end[i]=1;
 539     }
 540   }
 541
 542
 543
 544   if(firstLay==-1 || (secondLay==-1 && !onePoint)) return 0;
 545
 546   TClonesArray &arrSA= *fListOfUTracks;
 547   Int_t nFoundTracks=0;
 548
 549
 550   for(Int_t l0=0;l0<end[0];l0++){ //loop on layer 1
 551     AliITSRecPoint* cl0 = (AliITSRecPoint*)listlayer[0][l0];
 552     for(Int_t l1=0;l1<end[1];l1++){ //loop on layer 2
 553       AliITSRecPoint* cl1 = (AliITSRecPoint*)listlayer[1][l1];
 554       for(Int_t l2=0;l2<end[2];l2++){  //loop on layer 3
 555         AliITSRecPoint* cl2 = (AliITSRecPoint*)listlayer[2][l2];
 556         for(Int_t l3=0;l3<end[3];l3++){ //loop on layer 4
 557           AliITSRecPoint* cl3 = (AliITSRecPoint*)listlayer[3][l3];
 558           for(Int_t l4=0;l4<end[4];l4++){ //loop on layer 5
 559             AliITSRecPoint* cl4 = (AliITSRecPoint*)listlayer[4][l4];
 560             for(Int_t l5=0;l5<end[5];l5++){ //loop on layer 6
 561               AliITSRecPoint* cl5 = (AliITSRecPoint*)listlayer[5][l5];
 562               for(Int_t l6=0;l6<end[6];l6++){//loop on layer 7
 563                 AliITSRecPoint* cl6 = (AliITSRecPoint*)listlayer[6][l6];
 564                 for(Int_t l7=0;l7<end[7];l7++){ //loop on layer 8
 565                   AliITSRecPoint* cl7 = (AliITSRecPoint*)listlayer[7][l7];
 566                   Double_t x1,y1,z1,sx1,sy1,sz1;
 567                   Double_t x2,y2,z2,sx2,sy2,sz2;
 568                   AliITSRecPoint* p1=0;
 569                   AliITSRecPoint* p2=0;
 570                   Int_t index1=0,index2=0;
 571                   Int_t mrk1=0,mrk2=0;
 572                   switch(firstLay) {
 573                   case 0:
 574                     p1=cl0;
 575                     index1=clind[0][l0];mrk1=clmark[0][l0];
 576                     break;
 577                   case 1:
 578                     p1=cl1;
 579                     index1=clind[1][l1];mrk1=clmark[1][l1];
 580                     break;
 581                   case 2:
 582                     p1=cl2;
 583                     index1=clind[2][l2];mrk1=clmark[2][l2];
 584                     break;
 585                   case 3:
 586                     p1=cl3;
 587                     index1=clind[3][l3];mrk1=clmark[3][l3];
 588                     break;
 589                   case 4:
 590                     p1=cl4;
 591                     index1=clind[4][l4];mrk1=clmark[4][l4];
 592                     break;
 593                   case 5:
 594                     p1=cl5;
 595                     index1=clind[5][l5];mrk1=clmark[5][l5];
 596                     break;
 597                   case 6:
 598                     p1=cl6;
 599                     index1=clind[6][l6];mrk1=clmark[6][l6];
 600                     break;
 601                   }
 602                   switch(secondLay) {
 603                   case 1:
 604                     p2=cl1;
 605                     index2=clind[1][l1];mrk2=clmark[1][l1];
 606                     break;
 607                   case 2:
 608                     p2=cl2;
 609                     index2=clind[2][l2];mrk2=clmark[2][l2];
 610                     break;
 611                   case 3:
 612                     p2=cl3;
 613                     index2=clind[3][l3];mrk2=clmark[3][l3];
 614                     break;
 615                   case 4:
 616                     p2=cl4;
 617                     index2=clind[4][l4];mrk2=clmark[4][l4];
 618                     break;
 619                   case 5:
 620                     p2=cl5;
 621                     index2=clind[5][l5];mrk2=clmark[5][l5];
 622                     break;
 623                   case 6:
 624                     p2=cl6;
 625                     index2=clind[6][l6];mrk2=clmark[6][l6];
 626                     break;
 627                   case 7:
 628                     p2=cl7;
 629                     index2=clind[7][l7];mrk2=clmark[7][l7];
 630                     break;
 631
 632                   default:
 633                     p2=0;
 634                     index2=-1;mrk2=-1;
 635                     break;
 636                   }
 637                   Int_t layer = 0;
 638                   Double_t radius = 0.;
 639                   Double_t yclu1 = p1->GetY();
 640                   Double_t zclu1 = p1->GetZ();
 641                   layer=p1->GetLayer();
 642                   radius = fSegmentation->GetRadius(layer);
 643                   Double_t cv=0.,tgl2=0.,phi2=0.;
 644                   Int_t cln1=mrk1;
 645                   AliITSclusterTable* arr1 = (AliITSclusterTable*)GetClusterCoord(firstLay,cln1);
 646                   x1 = arr1->GetX();
 647                   y1 = arr1->GetY();
 648                   z1 = arr1->GetZ();
 649                   sx1 = arr1->GetSx();
 650                   sy1 = arr1->GetSy();
 651                   sz1 = arr1->GetSz();
 652
 653                   if(secondLay>0) {
 654                     Int_t cln2=mrk2;
 655                     AliITSclusterTable* arr2 = (AliITSclusterTable*)GetClusterCoord(secondLay,cln2);
 656                     x2 = arr2->GetX();
 657                     y2 = arr2->GetY();
 658                     z2 = arr2->GetZ();
 659                     sx2 = arr2->GetSx();
 660                     sy2 = arr2->GetSy();
 661                     sz2 = arr2->GetSz();
 662                     cv = Curvature(primaryVertex[0],primaryVertex[1],x1,y1,x2,y2);
 663                     tgl2 = (z2-z1)/TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
 664                     phi2 = TMath::ATan2((y2-y1),(x2-x1));
 665                   }
 666                   Double_t x,y,z = 0.;
 667                   Double_t xz[2];
 668                   xz[0]= p1->GetDetLocalX();
 669                   xz[1]= p1->GetDetLocalZ();
 670                   Bool_t check2;
 671                   check2 = fSegmentation->DetToGlobal(layer,xz[0], xz[1],x,y,z);
 672
 673
 674                   Double_t phiclrad;
 675                   phiclrad=TMath::ATan2(y,x);
 676                   if (phiclrad<0) phiclrad+=TMath::TwoPi();
 677                   else if (phiclrad>=TMath::TwoPi()) phiclrad-=TMath::TwoPi();
 678
 679                    Double_t alpha = fSegmentation->GetAlpha(fSegmentation->GetModule(phiclrad));
 680
 681                   AliITStrackU trac(alpha,radius,yclu1,zclu1,phi2,tgl2,cv,1,fNLayers);
 682                   if(cl7!=0){
 683                     trac.AddClusterV2(7,(clind[7][l7] & 0x0fffffff)>>0);
 684                     trac.AddClusterMark(7,clmark[7][l7]);
 685                   }
 686                   if(cl6!=0){
 687                     trac.AddClusterV2(6,(clind[6][l6] & 0x0fffffff)>>0);
 688                     trac.AddClusterMark(6,clmark[6][l6]);
 689                   }
 690                   if(cl5!=0) {
 691                     trac.AddClusterV2(5,(clind[5][l5] & 0x0fffffff)>>0);
 692                     trac.AddClusterMark(5,clmark[5][l5]);
 693                   }
 694                   if(cl4!=0){
 695                     trac.AddClusterV2(4,(clind[4][l4] & 0x0fffffff)>>0);
 696                     trac.AddClusterMark(4,clmark[4][l4]);
 697                   }
 698                   if(cl3!=0){
 699                     trac.AddClusterV2(3,(clind[3][l3] & 0x0fffffff)>>0);
 700                     trac.AddClusterMark(3,clmark[3][l3]);
 701                   }
 702                   if(cl2!=0){
 703                     trac.AddClusterV2(2,(clind[2][l2] & 0x0fffffff)>>0);
 704                     trac.AddClusterMark(2,clmark[2][l2]);
 705                   }
 706                   if(cl1!=0){
 707                     trac.AddClusterV2(1,(clind[1][l1] & 0x0fffffff)>>0);
 708                     trac.AddClusterMark(1,clmark[1][l1]);
 709                   }
 710                   if(cl0!=0){
 711                     trac.AddClusterV2(0,(clind[0][l0] & 0x0fffffff)>>0);
 712                     trac.AddClusterMark(0,clmark[0][l0]);
 713                   }
 714
 715                   //fit with Kalman filter using RefitAtBase()
 716                   AliITStrackU ot(trac, kFALSE);
 717                   ot.ResetCovariance(10.);
 718                   ot.ResetClusters();
 719
 720                   // Propagate inside the innermost layer with a cluster
 721                   if(ot.Propagate(ot.GetX()-0.1*ot.GetX())) {
 722                     if(RefitAtBase(AliITSRecoParam::GetrInsideITSscreen(),&ot,inx)){ //fit from layer 1 to layer 6
 723                       AliITStrackU otrack2(ot,kFALSE);
 724                       otrack2.ResetCovariance(10.);
 725                       otrack2.ResetClusters();
 726
 727                       //fit from last layer to layer 1
 728                       if(RefitAtBase(/*AliITSRecoParam::GetrInsideSPD1()*/fSegmentation->GetRadius(0)-0.3,&otrack2,inx)) {//check clind
 729                         new(arrSA[nFoundTracks]) AliITStrackU(otrack2);
 730                         ++nFoundTracks;
 731                       }
 732
 733                     }
 734                   }
 735
 736
 737                 }//end loop layer 8
 738               }//end loop layer 7
 739             }//end loop layer 6
 740           }//end loop layer 5
 741         }//end loop layer 4
 742       }//end loop layer 3
 743     }//end loop layer 2
 744   }//end loop layer 1
 745
 746
 747
 748
 749   if(fListOfUTracks->GetEntries()==0) return 0;
 750
 751   Int_t lowchi2 = FindTrackLowChiSquare();
 752   AliITStrackU* otrack = (AliITStrackU*)fListOfUTracks->At(lowchi2);
 753
 754   if(otrack==0) {
 755     return 0;
 756   }
 757
 758   CookLabel(otrack,0.); //MI change - to see fake ratio
 759   Int_t label=FindLabel(otrack);
 760   otrack->SetLabel(label);
 761
 762   Int_t indexc[fNLayers];
 763   for(Int_t i=0;i<fNLayers;i++) indexc[i]=0;
 764   for(Int_t nind=0;nind<otrack->GetNumberOfClusters();nind++){
 765     indexc[nind] = otrack->GetClusterIndex(nind);
 766   }
 767   //remove clusters of found track
 768   for(Int_t nlay=0;nlay<fNLayers;nlay++){
 769     for(Int_t cln=0;cln<otrack->GetNumberOfMarked(nlay);cln++){
 770       Int_t index = otrack->GetClusterMark(nlay,cln);
 771       fCluLayer[nlay]->RemoveAt(index);
 772       RemoveClusterCoord(nlay,index);
 773       fCluLayer[nlay]->Compress();
 774     }
 775   }
 776   return otrack;
 777
 778 }
 779
 780 //_______________________________________________________
 781 void AliITStrackerUpgrade::StoreTrack(AliITStrackV2 *t,AliESDEvent *event, Bool_t pureSA) /*const*/
 782 {
 783   //
 784   // Add new track to the ESD
 785   //
 786   AliESDtrack outtrack;
 787   outtrack.UpdateTrackParams(t,AliESDtrack::kITSin);
 788   if(pureSA) outtrack.SetStatus(AliESDtrack::kITSpureSA);
 789   outtrack.SetStatus(AliESDtrack::kITSout);
 790   outtrack.SetStatus(AliESDtrack::kITSrefit);
 791   Double_t sdedx[4]={0.,0.,0.,0.};
 792   for(Int_t i=0; i<4; i++) sdedx[i]=t->GetSampledEdx(i);
 793   outtrack.SetITSdEdxSamples(sdedx);
 794   event->AddTrack(&outtrack);
 795   return;
 796 }
 797
 798
 799 //_______________________________________________________
 800 Int_t AliITStrackerUpgrade::SearchClusters(Int_t layer,Double_t phiwindow,Double_t lambdawindow, AliITStrackU* trs,Double_t /*zvertex*/,Int_t pflag){
 801   //function used to to find the clusters associated to the track
 802
 803   AliDebug(2,"Starting...");
 804   if(ForceSkippingOfLayer(layer)) {
 805     AliDebug(2,Form("Forcing skipping of layer %i. Exiting",layer));
 806     return 0;
 807   }
 808
 809   Int_t nc=0;
 810   Double_t r=fSegmentation->GetRadius(layer);
 811   if(pflag==1){
 812     Float_t cx1,cx2,cy1,cy2;
 813     FindEquation(fPoint1[0],fPoint1[1],fPoint2[0],fPoint2[1],fPoint3[0],fPoint3[1],fCoef1,fCoef2,fCoef3);
 814     if (FindIntersection(fCoef1,fCoef2,fCoef3,-r*r,cx1,cy1,cx2,cy2)==0)
 815       return 0;
 816     Double_t fi1=TMath::ATan2(cy1-fPoint1[1],cx1-fPoint1[0]);
 817     Double_t fi2=TMath::ATan2(cy2-fPoint1[1],cx2-fPoint1[0]);
 818     fPhiEstimate=ChoosePoint(fi1,fi2,fPhic);
 819   }
 820
 821
 822   Int_t ncl = fCluLayer[layer]->GetEntries();
 823   AliDebug(2,Form(" Number of clusters %i in layer %i.",ncl,layer));
 824   for (Int_t index=0; index<ncl; index++) {
 825     AliITSRecPointU *c = (AliITSRecPointU*)fCluLayer[layer]->At(index);
 826     if (!c) continue;
 827     if (c->GetQ()<=0) continue;
 828     if(layer>1 && c->GetQ()<=fMinQ) continue;
 829
 830     AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(layer,index);
 831
 832     Double_t phi;
 833     Double_t xT,yT;
 834     xT=arr->GetX();
 835     yT=arr->GetY();
 836     phi=arr->GetPhi();
 837     if (TMath::Abs(phi-fPhiEstimate)>phiwindow)
 838       {
 839         continue;
 840       }
 841
 842     Double_t lambda = arr->GetLambda();
 843     if (TMath::Abs(lambda-fLambdac)>lambdawindow) continue;
 844
 845     if(trs->GetNumberOfClustersU()==trs->GetMaxNumberOfClusters()) return 0;
 846     if(trs->GetNumberOfMarked(layer)==trs->GetMaxNMarkedPerLayer()) return 0;
 847     Int_t orind = arr->GetOrInd();
 848     trs->AddClusterU(layer,orind);
 849     trs->AddClusterMark(layer,index);
 850
 851     nc++;
 852     fLambdac=lambda;
 853     fPhiEstimate=phi;
 854
 855     fPointc[0]=arr->GetX();
 856     fPointc[1]=arr->GetY();
 857   }
 858   return nc;
 859 }
 860
 861 //________________________________________________________________
 862 Bool_t AliITStrackerUpgrade::SetFirstPoint(Int_t lay, Int_t clu, Double_t* primaryVertex){
 863   // Sets the first point (seed) for tracking
 864
 865   AliITSRecPointU* cl = (AliITSRecPointU*)fCluLayer[lay]->At(clu);
 866   if(!cl) return kFALSE;
 867   if (cl->GetQ()<=0) return kFALSE;
 868   if(lay>1 && cl->GetQ()<=fMinQ) return kFALSE;
 869   AliITSclusterTable* arr = (AliITSclusterTable*)GetClusterCoord(lay,clu);
 870   fPhic = arr->GetPhi();
 871   fLambdac = arr->GetLambda();
 872   fPhiEstimate = fPhic;
 873   fPoint1[0]=primaryVertex[0];
 874   fPoint1[1]=primaryVertex[1];
 875   fPoint2[0]=arr->GetX();
 876   fPoint2[1]=arr->GetY();
 877   return kTRUE;
 878 }
 879
 880 //________________________________________________________________
 881 void AliITStrackerUpgrade::UpdatePoints(){
 882   //update of points for the estimation of the curvature
 883
 884   fPoint2[0]=fPoint3[0];
 885   fPoint2[1]=fPoint3[1];
 886   fPoint3[0]=fPointc[0];
 887   fPoint3[1]=fPointc[1];
 888
 889
 890 }
 891
 892 //___________________________________________________________________
 893 Int_t AliITStrackerUpgrade::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){
 894
 895   //given (x,y) of three recpoints (in global coordinates)
 896   //returns the parameters a,b,c of circonference x*x + y*y +a*x + b*y +c
 897
 898   Float_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 899   if(den==0) return 0;
 900   a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 901   b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 902   c = -x1*x1-y1*y1-a*x1-b*y1;
 903   return 1;
 904 }
 905 //__________________________________________________________________________
 906 Int_t AliITStrackerUpgrade::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){
 907
 908   //Finds the intersection between the circonference of the track and the circonference centered in (0,0) represented by one layer
 909   //c2 is -rlayer*rlayer
 910
 911   if(a1==0) return 0;
 912   Double_t m = c2-c1;
 913   Double_t aA = (b1*b1)/(a1*a1)+1;
 914   Double_t bB = (-2*m*b1/(a1*a1));
 915   Double_t cC = c2+(m*m)/(a1*a1);
 916   Double_t dD = bB*bB-4*aA*cC;
 917   if(dD<0) return 0;
 918
 919   y1 = (-bB+TMath::Sqrt(dD))/(2*aA);
 920   y2 = (-bB-TMath::Sqrt(dD))/(2*aA);
 921   x1 = (c2-c1-b1*y1)/a1;
 922   x2 = (c2-c1-b1*y2)/a1;
 923
 924   return 1;
 925 }
 926 //____________________________________________________________________
 927 Double_t AliITStrackerUpgrade::Curvature(Double_t x1,Double_t y1,Double_t
 928                                          x2,Double_t y2,Double_t x3,Double_t y3){
 929
 930   //calculates the curvature of track
 931   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 932   if(den==0) return 0;
 933   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 934   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 935   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
 936   Double_t xc=-a/2.;
 937
 938   if((a*a+b*b-4*c)<0) return 0;
 939   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
 940   if(rad==0) return 0;
 941
 942   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
 943   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
 944   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
 945   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
 946
 947   return 1/rad;
 948
 949 }
 950
 951
 952 //____________________________________________________________________
 953 Double_t AliITStrackerUpgrade::ChoosePoint(Double_t p1, Double_t p2, Double_t pp){
 954
 955   //Returns the point closest to pp
 956
 957   Double_t diff1 = p1-pp;
 958   Double_t diff2 = p2-pp;
 959
 960   if(TMath::Abs(diff1)<TMath::Abs(diff2)) fPhiEstimate=p1;
 961   else fPhiEstimate=p2;
 962   return fPhiEstimate;
 963
 964 }
 965
 966
 967 //_________________________________________________________________
 968 Int_t AliITStrackerUpgrade::FindTrackLowChiSquare() const {
 969   // returns track with lowest chi square
 970   Int_t dim=fListOfUTracks->GetEntries();
 971   if(dim<=1) return 0;
 972   AliITStrackV2* trk = (AliITStrackV2*)fListOfUTracks->At(0);
 973   Double_t minChi2=trk->GetChi2();
 974   Int_t index=0;
 975   for(Int_t i=1;i<dim;i++){
 976     trk = (AliITStrackV2*)fListOfUTracks->At(i);
 977     Double_t chi2=trk->GetChi2();
 978     if(chi2<minChi2){
 979       minChi2=chi2;
 980       index=i;
 981     }
 982   }
 983   return index;
 984 }
 985
 986 //__________________________________________________________
 987 Int_t AliITStrackerUpgrade::FindLabel(AliITStrackV2* track) const {
 988   //
 989
 990   Int_t labl[fNLayers][3];
 991   Int_t cnts[fNLayers][3];
 992   for(Int_t j=0;j<fNLayers;j++){
 993     for(Int_t k=0;k<3;k++){
 994       labl[j][k]=-2;
 995       cnts[j][k]=1;
 996     }
 997   }
 998   Int_t iNotLabel=0;
 999   for(Int_t i=0;i<track->GetNumberOfClusters(); i++) {
1000     Int_t indexc = track->GetClusterIndex(i);
1001     AliITSRecPointU* cl = (AliITSRecPointU*)GetCluster(indexc);
1002     AliDebug(2,Form(" cluster index %i  ",indexc));
1003     Int_t iLayer=cl->GetLayer();
1004     for(Int_t k=0;k<3;k++){
1005       labl[iLayer][k]=cl->GetLabel(k);
1006       if(labl[iLayer][k]<0) iNotLabel++;
1007     }
1008   }
1009   if(iNotLabel==3*track->GetNumberOfClusters()) return -2;
1010
1011   for(Int_t j1=0;j1<fNLayers; j1++) {
1012     for(Int_t j2=0; j2<j1;  j2++){
1013       for(Int_t k1=0; k1<3; k1++){
1014         for(Int_t k2=0; k2<3; k2++){
1015           if(labl[j1][k1]>=0 && labl[j1][k1]==labl[j2][k2] && cnts[j2][k2]>0){
1016             cnts[j2][k2]++;
1017             cnts[j1][k1]=0;
1018           }
1019         }
1020       }
1021     }
1022   }
1023
1024
1025   Int_t cntMax=0;
1026   Int_t label=-1;
1027   for(Int_t j=0;j<fNLayers;j++){
1028     for(Int_t k=0;k<3;k++){
1029       if(cnts[j][k]>cntMax && labl[j][k]>=0){
1030         cntMax=cnts[j][k];
1031         label=labl[j][k];
1032       }
1033     }
1034   }
1035
1036   Int_t lflag=0;
1037   for(Int_t i=0;i<fNLayers;i++)
1038     if(labl[i][0]==label || labl[i][1]==label || labl[i][2]==label) lflag++;
1039   if(lflag<track->GetNumberOfClusters()) label = -label;
1040   return label;
1041 }
1042 //_____________________________________________________________________________
1043 void AliITStrackerUpgrade::SetCalculatedWindowSizes(Int_t n, Float_t phimin, Float_t phimax, Float_t lambdamin, Float_t lambdamax){
1044   // Set sizes of the phi and lambda windows used for track finding
1045   fNloop = n;
1046   if(fPhiWin) delete [] fPhiWin;
1047   if(fLambdaWin) delete [] fLambdaWin;
1048   fPhiWin = new Double_t[fNloop];
1049   fLambdaWin = new Double_t[fNloop];
1050   Float_t stepPhi=(phimax-phimin)/(Float_t)(fNloop-1);
1051   Float_t stepLambda=(lambdamax-lambdamin)/(Float_t)(fNloop-1);
1052   for(Int_t k=0;k<fNloop;k++){
1053     Float_t phi=phimin+k*stepPhi;
1054     Float_t lam=lambdamin+k*stepLambda;
1055     fPhiWin[k]=phi;
1056     fLambdaWin[k]=lam;
1057   }
1058 }
1059 //_____________________________________________________________________________
1060 void AliITStrackerUpgrade::SetFixedWindowSizes(Int_t n, Double_t *phi, Double_t *lam){
1061   // Set sizes of the phi and lambda windows used for track finding
1062   fNloop = n;
1063   if(phi){ // user defined values
1064     fPhiWin = new Double_t[fNloop];
1065     fLambdaWin = new Double_t[fNloop];
1066     for(Int_t k=0;k<fNloop;k++){
1067       fPhiWin[k]=phi[k];
1068       fLambdaWin[k]=lam[k];
1069     }
1070   }
1071   else {  // default values
1072
1073     Double_t phid[33]   = {0.002,0.003,0.004,0.0045,0.0047,
1074                            0.005,0.0053,0.0055,
1075                            0.006,0.0063,0.0065,0.007,0.0073,0.0075,0.0077,
1076                            0.008,0.0083,0.0085,0.0087,0.009,0.0095,0.0097,
1077                            0.01,0.0105,0.011,0.0115,0.012,0.0125,0.013,0.0135,0.0140,0.0145};
1078     Double_t lambdad[33] = {0.003,0.004,0.005,0.005,0.005,
1079                             0.005,0.005,0.006,
1080                             0.006,0.006,0.006,0.007,0.007,0.007,0.007,
1081                             0.007,0.007,0.007,0.007,0.007,0.007,0.007,
1082                             0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008,0.008};
1083
1084     if(fNloop!=33){
1085       fNloop = 33;
1086     }
1087
1088
1089     fPhiWin = new Double_t[fNloop];
1090     fLambdaWin = new Double_t[fNloop];
1091
1092     Double_t factor=AliITSReconstructor::GetRecoParam()->GetFactorSAWindowSizes();
1093
1094     for(Int_t k=0;k<fNloop;k++){
1095       fPhiWin[k]=phid[k]*factor;
1096       fLambdaWin[k]=lambdad[k]*factor;
1097     }
1098
1099   }
1100
1101 }
1102 //_______________________________________________________________________
1103 void AliITStrackerUpgrade::GetCoorAngles(AliITSRecPointU* cl,Double_t &phi,Double_t &lambda, Double_t &x, Double_t &y,Double_t &z,Double_t* vertex){
1104   //Returns values of phi (azimuthal) and lambda angles for a given cluster
1105
1106   Double_t xz[2];
1107   xz[0]= cl->GetDetLocalX();
1108   xz[1]= cl->GetDetLocalZ() ;
1109   Bool_t check2;
1110   Int_t ilayer;
1111   ilayer = cl->GetLayer();
1112   check2 = fSegmentation->DetToGlobal(ilayer,xz[0], xz[1],x,y,z);
1113
1114   if(x!=0 && y!=0)
1115     phi=TMath::ATan2(y-vertex[1],x-vertex[0]);
1116   lambda=TMath::ATan2(z-vertex[2],TMath::Sqrt((x-vertex[0])*(x-vertex[0])+(y-vertex[1])*(y-vertex[1])));
1117 }
1118
1119 //________________________________________________________________________
1120 void AliITStrackerUpgrade::GetCoorErrors(AliITSRecPointU* cl,Float_t &sx,Float_t &sy, Float_t &sz){
1121
1122   //returns sigmax, y, z of cluster in global coordinates
1123   /*
1124     Double_t rot[9];     fGeom->GetRotMatrix(module,rot);
1125     Int_t lay,lad,det;
1126     AliITSgeomTGeo::GetModuleId(module,lay,lad,det);
1127
1128     Double_t alpha=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
1129     Double_t phi=TMath::Pi()/2+alpha;
1130     if (lay==1) phi+=TMath::Pi();
1131
1132     Float_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
1133   */
1134   Float_t covm[6];
1135   cl->GetGlobalCov(covm);
1136   sx=TMath::Sqrt(covm[0]);
1137   sy=TMath::Sqrt(covm[3]);
1138   sz=TMath::Sqrt(covm[5]);
1139
1140
1141
1142   //  Float_t covm[6];
1143   //  cl->GetGlobalCov(covm);
1144   //  sx=12*1e-04;//TMath::Sqrt(covm[0]);
1145   //  sy=12*1e-04;//TMath::Sqrt(covm[3]);
1146   //  sz=40*1e-04;//TMath::Sqrt(covm[5]);
1147   /*
1148     sx = TMath::Sqrt(sp*sp*cl->GetSigmaY2());
1149     sy = TMath::Sqrt(cp*cp*cl->GetSigmaY2());
1150     sz = TMath::Sqrt(cl->GetSigmaZ2());
1151   */
1152 }
1153 //_____________________________________________________________
1154 void AliITStrackerUpgrade::UnloadClusters() {
1155   //--------------------------------------------------------------------
1156   //This function unloads ITS clusters
1157   //--------------------------------------------------------------------
1158   for (Int_t i=0; i<fNLayers; i++) fLayers[i]->ResetClusters();
1159 }
1160 //______________________________________________________________________
1161 Bool_t AliITStrackerUpgrade::RefitAtBase(Double_t xx,AliITStrackU *track,
1162                                          const Int_t *clusters)
1163 {
1164   //--------------------------------------------------------------------
1165   // Simplified version for ITS upgrade studies -- does not use module info
1166   //--------------------------------------------------------------------
1167   Int_t index[fNLayers];
1168   Int_t k;
1169   for (k=0; k<fNLayers; k++) index[k]=-1;
1170   //
1171   for (k=0; k<fNLayers; k++) {
1172     index[k]=clusters[k];
1173   }
1174
1175   ULong_t trStatus=0;
1176   if(track->GetESDtrack()) trStatus=track->GetStatus();
1177   Int_t innermostlayer=0;
1178   if(trStatus&AliESDtrack::kTPCin)  {
1179     innermostlayer=fNLayers-1;
1180     Double_t drphi = TMath::Abs(track->GetD(0.,0.));
1181     for(innermostlayer=0; innermostlayer<fNLayers; innermostlayer++) {
1182       if( (drphi < (fSegmentation->GetRadius(innermostlayer)+1.)) ||
1183           index[innermostlayer] >= 0 ) break;
1184     }
1185     AliDebug(2,Form(" drphi  %f  innermost %d",drphi,innermostlayer));
1186   }
1187
1188   Int_t from, to, step;
1189   if (xx > track->GetX()) {
1190     from=innermostlayer; to=fNLayers;
1191     step=+1;
1192   } else {
1193     from=fNLayers-1; to=innermostlayer-1;
1194     step=-1;
1195   }
1196   TString dir = (step>0 ? "outward" : "inward");
1197
1198   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
1199     Double_t r=0.;
1200     r=fSegmentation->GetRadius(ilayer);
1201
1202     if (step<0 && xx>r){
1203       break;
1204     }
1205
1206
1207     Double_t oldGlobXYZ[3];
1208
1209     if (!track->GetXYZ(oldGlobXYZ)) {
1210       return kFALSE;
1211     }
1212     // continue if we are already beyond this layer
1213     Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
1214     if(step>0 && oldGlobR > r){
1215       continue; // going outward
1216     }
1217     if(step<0 && oldGlobR < r) {
1218       continue; // going inward
1219     }
1220     Double_t phi,z;
1221     if (!track->GetPhiZat(r,phi,z)){
1222       return kFALSE;
1223     }
1224     // only for ITS-SA tracks refit
1225     track->SetCheckInvariant(kFALSE);
1226
1227     if (phi<0) phi+=TMath::TwoPi();//from 0 to 360 (rad)
1228     else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();//
1229
1230     Double_t alpha = fSegmentation->GetAlpha(fSegmentation->GetModule(phi));
1231
1232     if (!track->Propagate(alpha,r)) {
1233       return kFALSE;
1234     }
1235
1236     const AliITSRecPointU *clAcc=0;
1237     Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
1238
1239     Int_t idx=index[ilayer];
1240     if (idx>=0) { // cluster in this layer
1241       Int_t cli = idx&0x0fffffff;
1242       const AliITSRecPointU *cl=(AliITSRecPointU *)fLayers[ilayer]->GetCluster(cli);
1243
1244       if (cl) {
1245         Int_t cllayer = (idx & 0xf0000000) >> 28;;
1246         Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
1247         if (chi2<maxchi2) {
1248           clAcc=cl;
1249           maxchi2=chi2;
1250         } else {
1251           return kFALSE;
1252         }
1253       }
1254     }
1255
1256     if (clAcc) {
1257       if (!UpdateMI(track,clAcc,maxchi2,idx)){
1258         return kFALSE;
1259       }
1260       track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
1261     }
1262
1263
1264     // Correct for material of the current layer
1265     // cross material
1266     // add time if going outward
1267     if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)){
1268       return kFALSE;
1269     }
1270     track->SetCheckInvariant(kFALSE);
1271   } // end loop on layers
1272
1273   if (!track->PropagateTo(xx,0.,0.)){
1274     return kFALSE;
1275   }
1276   return kTRUE;
1277 }
1278
1279 //_____________________________________________________________________
1280 Int_t AliITStrackerUpgrade::UpdateMI(AliITStrackU* track, const AliITSRecPointU* cl,Double_t chi2,Int_t index) const
1281 {
1282   //
1283   // Update ITS track
1284   //
1285   Int_t layer = (index & 0xf0000000) >> 28;
1286   track->SetClIndex(layer, index);
1287   if (TMath::Abs(cl->GetQ())<1.e-13) return 0;  // ingore the "virtual" clusters
1288
1289
1290   // Take into account the mis-alignment (bring track to cluster plane)
1291   Double_t xTrOrig=track->GetX();
1292   Double_t trxyz[3]; track->GetXYZ(trxyz);
1293   Float_t xclu1 = cl->GetX();//upgrade clusters tracking coordinate
1294   //Float_t yclu1 = cl->GetY();
1295   //Float_t zclu1 = cl->GetZ();
1296   if (!track->Propagate(xTrOrig+xclu1)){
1297     return 0;
1298   }
1299   AliCluster c(*cl);
1300   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
1301   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
1302   c.SetSigmaYZ(track->GetSigmaYZ(layer));
1303
1304   // Bring the track back to detector plane in ideal geometry
1305   Int_t updated = track->UpdateMI(&c,chi2,index);
1306   if (!track->Propagate(xTrOrig)) {
1307     return 0;
1308   }
1309   if(!updated){
1310     AliDebug(2,"update failed");
1311   }
1312   return updated;
1313 }
1314 //____________________________________________________________________-
1315 Int_t AliITStrackerUpgrade::CorrectForLayerMaterial(AliITStrackU *t,
1316                                                     Int_t layerindex,
1317                                                     Double_t oldGlobXYZ[3],
1318                                                     TString direction) {
1319   //-------------------------------------------------------------------
1320   // Propagate beyond layer and correct for material
1321   // (material budget in different ways according to fUseTGeo value)
1322   // Add time if going outward (PropagateTo or PropagateToTGeo)
1323   //-------------------------------------------------------------------
1324
1325   // Define budget mode:
1326   // 0: material from AliITSRecoParam (hard coded)
1327   // 1: material from TGeo in stepsof X cm (on the fly)
1328   //    X = AliITSRecoParam::GetStepSizeTGeo()
1329   // 2: material from lut
1330   // 3: material from TGeo in one step (same for all hypotheses)
1331   Int_t mode;
1332   switch(fUseTGeo) {
1333   case 0:
1334     mode=0;
1335     break;
1336   case 1:
1337     mode=1;
1338     break;
1339   case 2:
1340     mode=2;
1341     break;
1342   case 3:
1343     if(fTrackingPhase.Contains("Clusters2Tracks"))
1344       { mode=3; } else { mode=1; }
1345     break;
1346   case 4:
1347     if(fTrackingPhase.Contains("Clusters2Tracks"))
1348       { mode=3; } else { mode=2; }
1349     break;
1350   default:
1351     mode=0;
1352     break;
1353   }
1354   if(fTrackingPhase.Contains("Default")) {
1355     mode=0;
1356   }
1357   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1358   Double_t r = 0.;
1359   r=fSegmentation->GetRadius(layerindex);
1360   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
1361   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
1362   Double_t xToGo;
1363   if (!t->GetLocalXat(rToGo,xToGo)) {
1364     return 0;
1365   }
1366   Int_t index=fNLayers*fCurrentEsdTrack+layerindex;
1367
1368
1369   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
1370   Int_t nsteps=1;
1371   // back before material (no correction)
1372   Double_t rOld,xOld;
1373   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
1374   if (!t->GetLocalXat(rOld,xOld)) return 0;
1375   if (!t->Propagate(xOld)) return 0;
1376
1377   switch(mode) {
1378   case 1:
1379     x0=21.82;
1380     xOverX0 = fSegmentation->GetThickness(layerindex)/x0;
1381
1382     lengthTimesMeanDensity = xOverX0*x0;
1383     lengthTimesMeanDensity *= dir;
1384     // Bring the track beyond the material
1385     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)){
1386       return 0;
1387     }
1388     break;
1389   case 0:
1390     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
1391     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1392     break;
1393   case 2:
1394     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
1395     xOverX0 = fxOverX0Layer[layerindex];
1396     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
1397     lengthTimesMeanDensity *= dir;
1398     // Bring the track beyond the material
1399     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1400     break;
1401   case 3:
1402     if(!fxOverX0LayerTrks || index<0 || index>=fNLayers*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
1403     if(fxOverX0LayerTrks[index]<0) {
1404       nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
1405       if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
1406       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1407                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
1408       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
1409       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
1410       return 1;
1411     }
1412     xOverX0 = fxOverX0LayerTrks[index];
1413     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
1414     lengthTimesMeanDensity *= dir;
1415     // Bring the track beyond the material
1416     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)){
1417       return 0;
1418     }
1419     break;
1420   }
1421
1422
1423   return 1;
1424 }
1425
1426
1427 //_____________________________________________________________________________
1428 Double_t AliITStrackerUpgrade::GetPredictedChi2MI(AliITStrackU* track, const AliITSRecPointU *cluster,Int_t layer)
1429 {
1430   //
1431   // Compute predicted chi2
1432   //
1433   Float_t erry,errz,covyz;
1434   Float_t theta = track->GetTgl();
1435   Float_t phi   = track->GetSnp();
1436   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
1437   GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
1438   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
1439   // Take into account the mis-alignment (bring track to cluster plane)
1440   Double_t xTrOrig=track->GetX();
1441   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
1442   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
1443   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
1444   // Bring the track back to detector plane in ideal geometry
1445   // [mis-alignment will be accounted for in UpdateMI()]
1446   if (!track->Propagate(xTrOrig)) return 1000.;
1447   Float_t ny=0.,nz=0.;
1448   //GetNTeor(layer,cluster,theta,phi,ny,nz);
1449   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
1450   if (delta>1){
1451     chi2+=0.5*TMath::Min(delta/2,2.);
1452     chi2+=2.*cluster->GetDeltaProbability();
1453   }
1454   //
1455   track->SetNy(layer,ny);
1456   track->SetNz(layer,nz);
1457   track->SetSigmaY(layer,erry);
1458   track->SetSigmaZ(layer, errz);
1459   track->SetSigmaYZ(layer,covyz);
1460
1461   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
1462   return chi2;
1463 }
1464 //________________________________________________________________
1465 Int_t AliITStrackerUpgrade::GetError(Int_t /*layer*/,
1466                                      const AliITSRecPointU *cl,
1467                                      Float_t /*tgl*/,Float_t /*tgphitr*/,Float_t /*expQ*/,
1468                                      Float_t &erry,Float_t &errz,Float_t &covyz,
1469                                      Bool_t /*addMisalErr*/)
1470 {
1471   //
1472   // Calculate cluster position error for the upgrade
1473   //
1474   Int_t retval=0;
1475   covyz=0.;
1476   retval = GetErrorOrigRecPoint(cl,erry,errz,covyz);
1477
1478
1479
1480   return retval;
1481
1482 }
1483 //____________________________________________________________________-
1484
1485 Int_t AliITStrackerUpgrade::GetErrorOrigRecPoint(const AliITSRecPointU *cl,
1486                                                  Float_t &erry,Float_t &errz,Float_t &covyz)
1487 {
1488   //
1489   // Calculate cluster position error (just take error from AliITSRecPoint)
1490   //
1491   erry   = TMath::Sqrt(cl->GetSigmaY2());
1492   errz   = TMath::Sqrt(cl->GetSigmaZ2());
1493   covyz  = cl->GetSigmaYZ();
1494   return 1;
1495 }
1496 //__________________________
1497 void AliITStrackerUpgrade::GetNTeor(Int_t layer,const AliITSRecPointU* /*cl*/,
1498                                     Float_t tgl,Float_t tgphitr,
1499                                     Float_t &ny,Float_t &nz)
1500 {
1501   //
1502   // Get "mean shape" (original parametrization from AliITStrackerMI)
1503   //
1504   tgl = TMath::Abs(tgl);
1505   tgphitr = TMath::Abs(tgphitr);
1506
1507   // SPD
1508   if (layer==0) {
1509     ny = 1.+tgphitr*3.2;
1510     nz = 1.+tgl*0.34;
1511     return;
1512   }
1513   if (layer==1) {
1514     ny = 1.+tgphitr*3.2;
1515     nz = 1.+tgl*0.28;
1516     return;
1517   }
1518   // SSD
1519   if (layer==4 || layer==5) {
1520     ny = 2.02+tgphitr*1.95;
1521     nz = 2.02+tgphitr*2.35;
1522     return;
1523   }
1524   // SDD
1525   ny  = 6.6-2.7*tgphitr;
1526   nz  = 2.8-3.11*tgphitr+0.45*tgl;
1527   return;
1528 }
1529 //_________________________________________________________________
1530 Int_t AliITStrackerUpgrade::PropagateBack(AliESDEvent *event) {
1531   //--------------------------------------------------------------------
1532   // This functions propagates reconstructed ITS tracks back
1533   // The clusters must be loaded !
1534   //--------------------------------------------------------------------
1535   fTrackingPhase="PropagateBack";
1536   Int_t nentr=event->GetNumberOfTracks();
1537   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
1538   Int_t ntrk=0;
1539   for (Int_t i=0; i<nentr; i++) {
1540     AliESDtrack *esd=event->GetTrack(i);
1541
1542     if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
1543     if (esd->GetStatus()&AliESDtrack::kITSout) continue;
1544
1545     /*
1546       AliITStrackMI *t=0;
1547       try {
1548       t=new AliITStrackMI(*esd);
1549       } catch (const Char_t *msg) {
1550       //Warning("PropagateBack",msg);
1551       delete t;
1552       continue;
1553       }
1554     */
1555     AliITStrackU *t = new AliITStrackU(*esd);
1556     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
1557
1558     ResetTrackToFollow(*t);
1559
1560     /*
1561     // propagate to vertex [SR, GSI 17.02.2003]
1562     // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
1563     if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
1564     if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
1565     fTrackToFollow.StartTimeIntegral();
1566     // from vertex to outside pipe
1567     CorrectForPipeMaterial(&fTrackToFollow,"outward");
1568     }*/
1569     // Start time integral and add distance from current position to vertex
1570     Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
1571     fTrackToFollow.GetXYZ(xyzTrk);
1572     Double_t dst2 = 0.;
1573     for (Int_t icoord=0; icoord<3; icoord++) {
1574       Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
1575       dst2 += di*di;
1576     }
1577     fTrackToFollow.StartTimeIntegral();
1578     fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
1579
1580     fTrackToFollow.ResetCovariance(10.);
1581     fTrackToFollow.ResetClusters();
1582     //
1583     Int_t inx[fNLayers];
1584     for (Int_t k=0; k<fNLayers; k++) inx[k]=-1;
1585     Int_t nclusters = t->GetNumberOfClusters();
1586     for(Int_t ncl=0;ncl<nclusters;ncl++){
1587       Int_t index = t-> GetClIndex(ncl);
1588       Int_t lay = (index & 0xf0000000) >> 28;
1589       inx[lay]=index;
1590     }
1591     //
1592     if (RefitAtBase(fSegmentation->GetRadius(fNLayers-1),&fTrackToFollow,inx)) {
1593       //fTrackToFollow.SetLabel(t->GetLabel());//
1594       //fTrackToFollow.CookdEdx();
1595       //CookLabel(&fTrackToFollow,0.); //For comparison only
1596       fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
1597       //UseClusters(&fTrackToFollow);
1598       ntrk++;
1599     }
1600     delete t;
1601   }
1602
1603   AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
1604
1605   fTrackingPhase="Default";
1606
1607   return 0;
1608 }
1609 //_________________________________________________________________________
1610 AliCluster *AliITStrackerUpgrade::GetCluster(Int_t index) const {
1611   //--------------------------------------------------------------------
1612   //       Return pointer to a given cluster
1613   //--------------------------------------------------------------------
1614   Int_t l=(index & 0xf0000000) >> 28;
1615   Int_t c=(index & 0x0fffffff) >> 0;
1616   AliDebug(2,Form("index %i  cluster index %i layer %i", index,c,l));
1617   return fLayers[l]->GetCluster(c);
1618 }
1619 //______________________________________________________________________________
1620 Int_t AliITStrackerUpgrade::CorrectForPipeMaterial(AliITStrackU *t, TString direction) {
1621   //-------------------------------------------------------------------
1622   // Propagate beyond beam pipe and correct for material
1623   // (material budget in different ways according to fUseTGeo value)
1624   // Add time if going outward (PropagateTo or PropagateToTGeo)
1625   //-------------------------------------------------------------------
1626
1627   // Define budget mode:
1628   // 0: material from AliITSRecoParam (hard coded)
1629   // 1: material from TGeo in one step (on the fly)
1630   // 2: material from lut
1631   // 3: material from TGeo in one step (same for all hypotheses)
1632   Int_t mode;
1633   switch(fUseTGeo) {
1634   case 0:
1635     mode=0;
1636     break;
1637   case 1:
1638     mode=1;
1639     break;
1640   case 2:
1641     mode=2;
1642     break;
1643   case 3:
1644     if(fTrackingPhase.Contains("Clusters2Tracks"))
1645       { mode=3; } else { mode=1; }
1646     break;
1647   case 4:
1648     if(fTrackingPhase.Contains("Clusters2Tracks"))
1649       { mode=3; } else { mode=2; }
1650     break;
1651   default:
1652     mode=0;
1653     break;
1654   }
1655   if(fTrackingPhase.Contains("Default")) mode=0;
1656   Int_t index=fCurrentEsdTrack;
1657
1658   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1659   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
1660   Double_t xToGo;
1661   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1662
1663   Double_t xOverX0,x0,lengthTimesMeanDensity;
1664
1665   switch(mode) {
1666   case 0:
1667     xOverX0 = AliITSRecoParam::GetdPipe();
1668     x0 = AliITSRecoParam::GetX0Be();
1669     lengthTimesMeanDensity = xOverX0*x0;
1670     lengthTimesMeanDensity *= dir;
1671     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1672     break;
1673   case 1:
1674     if (!t->PropagateToTGeo(xToGo,1)) return 0;
1675     break;
1676   case 2:
1677     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
1678     xOverX0 = fxOverX0Pipe;
1679     lengthTimesMeanDensity = fxTimesRhoPipe;
1680     lengthTimesMeanDensity *= dir;
1681     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1682     break;
1683   case 3:
1684     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
1685     if(fxOverX0PipeTrks[index]<0) {
1686       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1687       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1688                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
1689       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
1690       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
1691       return 1;
1692     }
1693     xOverX0 = fxOverX0PipeTrks[index];
1694     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
1695     lengthTimesMeanDensity *= dir;
1696     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1697     break;
1698   }
1699
1700   return 1;
1701 }
1702 //__________________________________________________________________________
1703 Int_t AliITStrackerUpgrade::RefitInward(AliESDEvent *event) {
1704   //--------------------------------------------------------------------
1705   // This functions refits ITS tracks using the
1706   // "inward propagated" TPC tracks
1707   // The clusters must be loaded !
1708   //--------------------------------------------------------------------
1709   fTrackingPhase="RefitInward";
1710   Int_t nentr=event->GetNumberOfTracks();
1711   for (Int_t i=0; i<nentr; i++) {
1712     AliESDtrack *esd=event->GetTrack(i);
1713
1714     // if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
1715     // if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
1716     //if (esd->GetStatus()&AliESDtrack::kTPCout)
1717     //  if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
1718     AliITStrackU *t=0;
1719     try {
1720       t=new AliITStrackU(*esd);
1721     } catch (const Char_t *msg) {
1722       //Warning("RefitInward",msg);
1723       delete t;
1724       continue;
1725     }
1726
1727     ResetTrackToFollow(*t);
1728     // fTrackToFollow.ResetClusters();
1729     //  if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
1730     //  fTrackToFollow.ResetCovariance(10.);
1731     //Refitting...
1732     //The beam pipe
1733     if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
1734       fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
1735       AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
1736       Double_t r[3]={0.,0.,0.};
1737       Double_t maxD=3.;
1738       esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
1739     }
1740     delete t;
1741   }
1742
1743   fTrackingPhase="Default";
1744   return 0;
1745 }
1746
1747
1748