ITS/UPGRADE/AliITSUTrackerSA.cxx

   1 //--------------------------------------------------------------------------------
   2 //               Implementation of the ITS tracker class
   3 //    It reads AliITSUClusterPix clusters and and fills the ESD with tracks
   4 //
   5 //    The algorithm implemented here takes inspiration from UniCA code of FIAS
   6 //    group.
   7 //--------------------------------------------------------------------------------
   8
   9 #include <TBranch.h>
  10 #include <TMath.h>
  11 using TMath::Abs;
  12 using TMath::Sort;
  13 using TMath::Sqrt;
  14 #include <TTree.h>
  15 #include <algorithm>
  16 using std::sort;
  17
  18 // Vc library
  19 //#include "Vc/Vc"
  20 //#include "AliITSUTrackerSAauxVc.h" // Structs and other stuff using Vc library
  21 #include "AliLog.h"
  22 #include "AliESDEvent.h"
  23 #include "AliITSUClusterPix.h"
  24 #include "AliITSUTrackerSA.h"
  25 #include "AliITSUReconstructor.h"
  26 #include "AliITSURecoDet.h"
  27 #include "AliESDtrack.h"
  28
  29 #include <Riostream.h>
  30
  31 using std::cout;
  32 using std::endl;
  33 using std::flush;
  34
  35 //#include "AliITSUtrackSA.h"      // Some dedicated SA track class ?
  36
  37 ClassImp(AliITSUTrackerSA)
  38
  39 const Double_t AliITSUTrackerSA::fgkToler =  1e-6;// tolerance for layer on-surface check
  40 const Double_t AliITSUTrackerSA::fgkChi2Cut =  600.f;
  41
  42 //________________________________________________________________________________
  43 AliITSUTrackerSA::AliITSUTrackerSA(AliITSUReconstructor* rec) :
  44 fReconstructor(rec),
  45 fITS(0),
  46 fMatLUT(0),
  47 fUseMatLUT(kFALSE),
  48 fCurrMass(0.14),
  49 //
  50 fClustersTC(),
  51 fChi2Cut( fgkChi2Cut ),
  52 fPhiCut( 1  ),
  53 fRPhiCut( 1 ),
  54 fZCut( 1 )
  55 {
  56   //--------------------------------------------------------------------
  57   // This default constructor needs to be provided
  58   //--------------------------------------------------------------------
  59   if (rec) Init(rec);
  60 }
  61
  62 //________________________________________________________________________________
  63 AliITSUTrackerSA::AliITSUTrackerSA(const AliITSUTrackerSA &t): AliTracker(t),
  64                                                                fReconstructor(t.fReconstructor),
  65                                                                fITS(t.fITS),
  66                                                                fMatLUT(t.fMatLUT),
  67                                                                fUseMatLUT(t.fUseMatLUT),
  68                                                                fCurrMass(t.fCurrMass),
  69   //
  70                                                                fClustersTC(),
  71                                                                fChi2Cut(fgkChi2Cut),
  72                                                                fPhiCut(),
  73                                                                fRPhiCut(),
  74                                                                fZCut()
  75 {
  76   //--------------------------------------------------------------------
  77   // The copy constructor is protected
  78   //--------------------------------------------------------------------
  79 }
  80
  81 //________________________________________________________________________________
  82 AliITSUTrackerSA::~AliITSUTrackerSA()
  83 {
  84   // d-tor
  85   delete fMatLUT;
  86 }
  87
  88
  89 //_________________________________________________________________________
  90 void AliITSUTrackerSA::Init(AliITSUReconstructor* rec)
  91 {
  92   // init with external reconstructor
  93   //
  94   fITS = rec->GetITSInterface();
  95   //
  96   // create material lookup table
  97   const int kNTest = 1000;
  98   const double kStepsPerCM=5;
  99   fMatLUT  = new AliITSUMatLUT(fITS->GetRMin(),fITS->GetRMax(),Nint(kStepsPerCM*(fITS->GetRMax()-fITS->GetRMin())));
 100   double zmn = 1e6;
 101   for (int ilr=fITS->GetNLayers();ilr--;) {
 102     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
 103     if (zmn>Abs(lr->GetZMin())) zmn = Abs(lr->GetZMin());
 104     if (zmn>Abs(lr->GetZMax())) zmn = Abs(lr->GetZMax());
 105   }
 106   fMatLUT->FillData(kNTest,-zmn,zmn);
 107   //
 108 }
 109
 110 //________________________________________________________________________________
 111 Int_t AliITSUTrackerSA::Clusters2Tracks(AliESDEvent *event) {
 112   //--------------------------------------------------------------------
 113   // This is the main tracking function
 114   // The clusters must already be loaded
 115   //--------------------------------------------------------------------
 116
 117   // Possibly, create the track "seeds" (combinatorial)
 118
 119   // Possibly, increment the seeds with additional clusters (Kalman)
 120
 121   // Possibly, (re)fit the found tracks
 122
 123   // Three iterations:
 124   // - High momentum first;
 125   // - Low momentum with vertex constraint;
 126   // - Everything else.
 127
 128   CellsCreation(0);
 129   CellularAutomaton(event);
 130   // VertexFinding();
 131   // CellsCreation(1);
 132   // CellularAutomaton(event);
 133
 134   return 0;
 135 }
 136
 137 //________________________________________________________________________________
 138 Int_t AliITSUTrackerSA::PropagateBack(AliESDEvent * event) {
 139   //--------------------------------------------------------------------
 140   // Here, we implement the Kalman smoother ?
 141   // The clusters must be already loaded
 142   //--------------------------------------------------------------------
 143   Int_t n=event->GetNumberOfTracks();
 144   Int_t ntrk=0;
 145   Int_t ngood=0;
 146   for (Int_t i=0; i<n; i++) {
 147     AliESDtrack *esdTrack=event->GetTrack(i);
 148
 149     if ((esdTrack->GetStatus()&AliESDtrack::kITSin)==0) continue;
 150
 151     AliITSUTrackCooked track(*esdTrack);
 152
 153     track.ResetCovariance(10.);
 154
 155     int points[2*AliITSUAux::kMaxLayers];
 156     for (UInt_t k=0; k<2*AliITSUAux::kMaxLayers; k++)
 157       points[k]=-1;
 158     Int_t nc=track.GetNumberOfClusters();
 159     for (Int_t k=0; k<nc; k++) {
 160       const int layer = (track.GetClusterIndex(k)&0xf0000000)>>28;
 161       const int idx = (track.GetClusterIndex(k)&0x0fffffff);
 162       points[layer<<1]=idx;
 163     }
 164
 165     if (RefitTrack(&track,points,40,1)>=0) {
 166
 167       CookLabel(&track, 0.); //For comparison only
 168       Int_t label=track.GetLabel();
 169       if (label>0) ngood++;
 170
 171       esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSout);
 172       ntrk++;
 173     }
 174   }
 175
 176   Info("PropagateBack","Back propagated tracks: %d",ntrk);
 177   if (ntrk)
 178     Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk);
 179
 180   return 0;
 181 }
 182
 183 //________________________________________________________________________________
 184 Int_t AliITSUTrackerSA::RefitInward(AliESDEvent * event) {
 185   //--------------------------------------------------------------------
 186   // Some final refit, after the outliers get removed by the smoother ?
 187   // The clusters must be loaded
 188   //--------------------------------------------------------------------
 189   Int_t n=event->GetNumberOfTracks();
 190   Int_t ntrk=0;
 191   Int_t ngood=0;
 192   for (Int_t i=0; i<n; i++) {
 193     AliESDtrack *esdTrack=event->GetTrack(i);
 194
 195     if ((esdTrack->GetStatus()&AliESDtrack::kITSout)==0) continue;
 196
 197     AliITSUTrackCooked track(*esdTrack);
 198
 199     track.ResetCovariance(10.);
 200
 201     int points[2*AliITSUAux::kMaxLayers];
 202     for (UInt_t k=0; k<2*AliITSUAux::kMaxLayers; k++)
 203       points[k]=-1;
 204     Int_t nc=track.GetNumberOfClusters();
 205     for (Int_t k=0; k<nc; k++) {
 206       const int layer = (track.GetClusterIndex(k)&0xf0000000)>>28;
 207       const int idx = (track.GetClusterIndex(k)&0x0fffffff);
 208       points[layer<<1]=idx;
 209     }
 210
 211     if (RefitTrack(&track,points,1.8,1)>=0) { //2.1,1)>=0) {
 212
 213       //if (!track.PropagateTo(1.8, 2.27e-3, 35.28*1.848)) continue;
 214       CookLabel(&track, 0.); //For comparison only
 215       Int_t label=track.GetLabel();
 216       if (label>0) ngood++;
 217
 218       //cout << esdTrack->GetStatus() << " ";
 219       esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSrefit);
 220       //cout << esdTrack->GetStatus() << endl;
 221       ntrk++;
 222     }
 223   }
 224
 225   Info("RefitInward","Refitted tracks: %d",ntrk);
 226   if (ntrk)
 227     Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood)/ntrk);
 228
 229   return 0;
 230 }
 231
 232 //________________________________________________________________________________
 233 Int_t AliITSUTrackerSA::LoadClusters(TTree *cluTree) {
 234   //--------------------------------------------------------------------
 235   // This function reads the ITSU clusters from the tree,
 236   // sort them, distribute over the internal tracker arrays, etc
 237   //--------------------------------------------------------------------
 238   fITS->LoadClusters(cluTree);
 239   fITS->ProcessClusters();
 240   //
 241   for(int iL=0; iL<7; ++iL) {
 242     fClustersTC[iL]=*fITS->GetLayerActive(iL)->GetClustersAddress();
 243     AliITSURecoLayer* lr = fITS->GetLayerActive(iL) ; // assign the layer which the cluster belongs to
 244     for(int iC=0;iC<fClustersTC[iL]->GetEntriesFast();++iC) {
 245       const AliITSUClusterPix *cl = (AliITSUClusterPix*)fClustersTC[iL]->At(iC);
 246       float pos[3];
 247       cl->GetGlobalXYZ(pos);
 248       float phi = TMath::PiOver2();
 249       if(Abs(pos[0])>1e-9) {
 250         phi=TMath::ATan2(pos[1]-GetY(),pos[0]-GetX());
 251         if(phi<0.f) phi+=TMath::TwoPi();
 252       } else if(pos[1]<0.f) phi *= 3.f ;
 253       AliITSURecoSens* sens = lr->GetSensorFromID(cl->GetVolumeId());
 254       const float alpha = sens->GetPhiTF();
 255       const float cov[3]={cl->GetSigmaZ2(),cl->GetSigmaYZ(),cl->GetSigmaY2()};
 256
 257       fLayer[iL].AddPoint(pos,cov,phi,alpha);
 258       for ( int i=0 ; i<3; ++i ) {
 259         fLayer[iL].Points.back().Label[i] = (cl->GetLabel(i)<0) ? -1 : cl->GetLabel(i);
 260       }
 261     }
 262
 263     fLayer[iL].Sort();
 264
 265   }
 266   return 0;
 267 }
 268
 269 //________________________________________________________________________________
 270 void AliITSUTrackerSA::UnloadClusters() {
 271   //--------------------------------------------------------------------
 272   // This function unloads ITSU clusters from the RAM
 273   //--------------------------------------------------------------------
 274   for(int i=0;i<7;++i)
 275     fLayer[i].Clear();
 276   for(int i=0;i<5;++i)
 277     fCells[i].clear();
 278
 279 }
 280
 281 //________________________________________________________________________________
 282 void AliITSUTrackerSA::CellularAutomaton(AliESDEvent *event) {
 283
 284   // Here it's implemented the Cellular Automaton routine
 285   // Firstly the level of each doublet is set according to the level of
 286   // the neighbour doublets.
 287   // Doublet are considered to be neighbour if they share one point and the
 288   // phi and theta direction difference of the two is below a cut value.
 289   Int_t ntrk=0,ngood=0;
 290
 291   for( int iL = 1; iL < 5; ++iL ) {
 292     for( size_t iC1 = 0; iC1 < fCells[iL-1].size(); ++iC1 ) {
 293       for( size_t iC2 = 0; iC2 < fCells[iL].size(); ++iC2 ) {
 294         if( fCells[iL-1][iC1].Points[1]==fCells[iL][iC2].Points[0] &&
 295             fCells[iL-1][iC1].Points[2]==fCells[iL][iC2].Points[1] &&
 296             fCells[iL-1][iC1].Level >= fCells[iL][iC2].Level - 1 ) {
 297           // The implementation of the curvature based matching has to be studied.
 298           fCells[iL][iC2].Level = fCells[iL-1][iC1].Level+1;
 299           fCells[iL][iC2].Neighbours.push_back(iC1);
 300         }
 301       }
 302     }
 303   }
 304
 305   for (int level = 5; level > 0; --level ) {
 306     vector<Road> roads;
 307     roads.reserve(100); // should reserve() be based on number of clusters on outermost layer?
 308     for (int iCL=4; iCL >= level-1; --iCL ) {
 309       for (size_t iCell = 0; iCell < fCells[iCL].size(); ++iCell) {
 310         if ( fCells[iCL][iCell].Level != level ) continue;
 311         roads.push_back( Road(iCL,iCell) );
 312         for( size_t iN=0;iN<fCells[iCL][iCell].Neighbours.size(); ++iN ) {
 313           const int currD = iCL - 1;
 314           const int neigh = fCells[iCL][iCell].Neighbours[iN];
 315           if( iN > 0 ) roads.push_back(roads.back());
 316           CandidatesTreeTraversal(roads,neigh,currD);
 317         }
 318         fCells[iCL][iCell].Level = -1;
 319       }
 320     }
 321
 322     // for(size_t iR=0; iR<roads.size(); ++iR) {
 323     //   cout << "ROAD " << iR << " | ";
 324     //   for(int i=0;i<5;++i) {
 325     //     if(roads[iR][i]<0) continue;
 326     //     else {
 327     //       if(roads[iR].Label==-1){
 328     //         roads[iR].Label = fCells[i][roads[iR][i]].Label;
 329     //         if(roads[iR].Label==-1) roads[iR].Label--;
 330     //       }
 331     //       if (fCells[i][roads[iR][i]].Label!=roads[iR].Label&&roads[iR].Label>-1) {
 332     //         roads[iR].Label = -1;
 333     //         if(fCells[i][roads[iR][i]].Label==-1) roads[iR].Label--;
 334     //       }
 335
 336     //       cout << fCells[i][roads[iR][i]].Label << " ";
 337     //     }
 338     //   }
 339     //   cout << " | " << roads[iR].Label << " | " << roads[iR].N << endl;
 340     // }
 341     vector<AliITSUTrackCooked> candidates;
 342     candidates.reserve(roads.size());
 343
 344     for (size_t iR=0; iR<roads.size(); ++iR) {
 345       if(roads[iR].N != level) {
 346         continue;
 347       }
 348
 349       int points[2*AliITSUAux::kMaxLayers];
 350       for(unsigned int i=0;i<2*AliITSUAux::kMaxLayers;++i) points[i] = -1;
 351       for(int i=0;i<5;++i) {
 352         if(roads[iR].Elements[i]<0) continue;
 353         points[( i )<<1]=fLayer[ i ](fCells[i][roads[iR].Elements[i]].Points[0]);
 354         points[(i+1)<<1]=fLayer[i+1](fCells[i][roads[iR].Elements[i]].Points[1]);
 355         points[(i+2)<<1]=fLayer[i+2](fCells[i][roads[iR].Elements[i]].Points[2]);
 356       }
 357
 358       candidates.push_back(AliITSUTrackCooked());
 359
 360       InitTrackParams(candidates.back(),points);
 361       const double chi2 = RefitTrack( (AliExternalTrackParam*)&candidates.back(), points, 0. ,-1);
 362
 363       if ( chi2 < 0. ) {
 364         // cout << "FAIL: " << chi2 << endl;
 365         // for(unsigned int i=0;i<2*AliITSUAux::kMaxLayers;++i)
 366         //   cout << points[i] << " ";
 367         // cout << endl;
 368         candidates.back().SetChi2( 1e27 );
 369       } else candidates.back().SetChi2( chi2 );
 370       candidates.back().SetLabel(roads[iR].Label);
 371     }
 372
 373     vector<int> index;
 374     index.reserve(candidates.size());
 375     for ( size_t i = 0; i < candidates.size(); ++i ) index.push_back(i);
 376     Comparison<AliITSUTrackCooked> comp(&candidates);
 377     sort(index.begin(),index.end(),comp);
 378
 379     for ( size_t cand = 0; cand < candidates.size(); ++cand ) {
 380       const int ii = index[cand];
 381
 382       if ( candidates[ii].GetChi2() < 0. ) continue;
 383
 384       // cout << candidates[ii].GetChi2() << " " << candidates[ii].GetNumberOfClusters() << " | " << candidates[ii].GetLabel() << " | ";
 385       // for(int i=0;i<candidates[ii].GetNumberOfClusters();++i) {
 386       //   cout<< GetCluster(candidates[ii].GetClusterIndex(i))->GetLabel(0) << " ";
 387       // }
 388       // cout << endl;
 389
 390       if( candidates[ii].GetChi2()/candidates[ii].GetNumberOfClusters() > fgkChi2Cut ) {
 391         break;
 392       }
 393       bool goodTrack = true;
 394       for ( int point = 0; point < candidates[ii].GetNumberOfClusters(); ++point ) {
 395         int layer = (candidates[ii].GetClusterIndex(point)&0xf0000000)>>28;
 396         int ind = (candidates[ii].GetClusterIndex(point)&0x0fffffff);
 397
 398         if( (fLayer[ layer ].Points[ ind ].Used ) ) {
 399           goodTrack = false;
 400         }
 401
 402       }
 403       if(!goodTrack) {
 404         continue;
 405       }
 406       for ( int point = 0; point < candidates[ii].GetNumberOfClusters(); ++point ) {
 407         int layer = (candidates[ii].GetClusterIndex(point)&0xf0000000)>>28;
 408         int ind = (candidates[ii].GetClusterIndex(point)&0x0fffffff);
 409         fLayer[ layer ].Points[ ind ].Used = true;
 410       }
 411
 412       AliESDtrack outTrack;
 413       CookLabel((AliKalmanTrack*)&candidates[ii],0.f);
 414       ntrk++;
 415       if(candidates[ii].GetChi2()>0) ngood++;
 416
 417       // cout << candidates[ii].GetChi2() << " " << candidates[ii].GetNumberOfClusters() << " | " << candidates[ii].GetLabel() << " | ";
 418       // for(int i=0;i<candidates[ii].GetNumberOfClusters();++i) {
 419       //   cout<< GetCluster(candidates[ii].GetClusterIndex(i))->GetLabel(0) << " ";
 420       // }
 421       // cout << endl;
 422
 423       outTrack.UpdateTrackParams((AliKalmanTrack*)&candidates[ii],AliESDtrack::kITSin);
 424       outTrack.SetLabel(candidates[ii].GetLabel());
 425       event->AddTrack(&outTrack);
 426     }
 427   }
 428   Info("Clusters2Tracks","Reconstructed tracks: %d",ntrk);
 429   if (ntrk)
 430     Info("Clusters2Tracks","Good tracks/reconstructed: %f",Float_t(ngood)/ntrk);
 431 }
 432
 433
 434 //________________________________________________________________________________
 435 void AliITSUTrackerSA::CellsCreation(const int &cutLevel) {
 436   // Make associations between two points on adjacent layers within an azimuthal window.
 437   // Under consideration:
 438   // - track parameter estimation using the primary vertex position
 439   // To do:
 440   // - last iteration
 441
 442   float phiCut = 7.f;
 443   if( cutLevel==0 ) phiCut = fPhiCut;
 444
 445   // Doublets creation
 446   vector<Cell> doublets[6];
 447   for( int iL = 0 ; iL < 6 ; ++iL ) {
 448     for ( int iC1 = 0 ; iC1 < fLayer[iL].N ; ++iC1 ) {
 449       for ( int iC2 = 0; iC2 < fLayer[iL+1].N ; ++iC2 ) {
 450         const float dPhi = Abs( fLayer[iL][iC1].Phi - fLayer[iL+1][iC2].Phi );
 451         if( dPhi < phiCut || Abs( dPhi - TMath::TwoPi() ) < phiCut) {
 452           doublets[iL].push_back(Cell(iC1,iC2));
 453           if(Abs(fLayer[iL][iC1].XYZ[0]-fLayer[iL+1][iC2].XYZ[0])<1e-32) {
 454             doublets[iL].back().Param[0] = 1e32;
 455           } else {
 456             doublets[iL].back().Param[0] = (fLayer[iL][iC1].XYZ[1]-fLayer[iL+1][iC2].XYZ[1])/(fLayer[iL][iC1].XYZ[0]-fLayer[iL+1][iC2].XYZ[0]);
 457           }
 458           const float r1  = Sqrt(fLayer[iL][iC1].XYZ[0] * fLayer[iL][iC1].XYZ[0] + fLayer[iL][iC1].XYZ[1] * fLayer[iL][iC1].XYZ[1]);
 459           const float r2  = Sqrt(fLayer[iL+1][iC2].XYZ[0] * fLayer[iL+1][iC2].XYZ[0] + fLayer[iL+1][iC2].XYZ[1] * fLayer[iL+1][iC2].XYZ[1]);
 460           doublets[iL].back().Param[1] = (fLayer[iL][iC1].XYZ[2]-fLayer[iL+1][iC2].XYZ[2])/(r1-r2);
 461           doublets[iL].back().Label=-1;
 462           for(int i=0;i<3;++i) {
 463             for(int j=0;j<3;++j) {
 464               if(fLayer[iL][iC1].Label[i]>-1&&fLayer[iL][iC1].Label[i]==fLayer[iL+1][iC2].Label[j])
 465                 doublets[iL].back().Label = fLayer[iL][iC1].Label[i];
 466             }
 467           }
 468         } else if( fLayer[iL+1][iC2].Phi - fLayer[iL][iC1].Phi > phiCut ) break;
 469       }
 470
 471     }
 472   }
 473
 474   // Triplets creation
 475   for( int iL = 5; iL > 0; --iL ) {
 476     fCells[iL-1].clear();
 477     for ( size_t iD2 = 0; iD2 < doublets[iL].size(); ++iD2 ) {
 478       for ( size_t iD1 = 0; iD1 < doublets[iL-1].size(); ++iD1 ) {
 479         const int id1 = doublets[iL-1][iD1].Points[1];
 480         const int id2 = doublets[iL][iD2].Points[0];
 481         if ( id1 == id2 ) {
 482           const int id3 = doublets[iL][iD2].Points[1];
 483           const float r3 = Sqrt( fLayer[iL+1][id3].XYZ[0] * fLayer[iL+1][id3].XYZ[0] + fLayer[iL+1][id3].XYZ[1] * fLayer[iL+1][id3].XYZ[1] );
 484           const float r2 = Sqrt( fLayer[iL][id2].XYZ[0] * fLayer[iL][id2].XYZ[0] + fLayer[iL][id2].XYZ[1] * fLayer[iL][id2].XYZ[1] );
 485           const float extrZ3 = doublets[iL-1][iD1].Param[1] * ( r3 - r2 ) + fLayer[iL][id2].XYZ[2] ;
 486           const int iii = doublets[iL-1][iD1].Points[0];
 487           if ( Abs ( extrZ3 - fLayer[iL+1][id3].XYZ[2] ) < fZCut ) {
 488             fCells[iL-1].push_back(Cell(doublets[iL-1][iD1].Points[0],id2,id3));
 489             fCells[iL-1].back().Param[0] = Curvature(fLayer[iL+1][id3].XYZ[0],fLayer[iL+1][id3].XYZ[1],fLayer[iL][id2].XYZ[0],fLayer[iL][id2].XYZ[1],fLayer[iL-1][iii].XYZ[0],fLayer[iL-1][iii].XYZ[1]);
 490             fCells[iL-1].back().Param[1] = doublets[iL][iD2].Param[1];
 491             if(doublets[iL-1][iD1].Label==doublets[iL][iD2].Label&&doublets[iL][iD2].Label!=-1)
 492               fCells[iL-1].back().Label=doublets[iL][iD2].Label;
 493             else
 494               fCells[iL-1].back().Label=-1;
 495           }
 496         }
 497       }
 498     }
 499   }
 500
 501 }
 502
 503 //______________________________________________________________________________
 504 Bool_t AliITSUTrackerSA::InitTrackParams(AliITSUTrackCooked &track, int points[])
 505 {
 506   // Set the initial guess on track kinematics for propagation.
 507   // Assume at least 3 points available
 508   int lrOcc[AliITSUAux::kMaxLayers], nCl=0;
 509   //
 510   // we will need endpoints and middle layer
 511   for (int i=fITS->GetNLayersActive()-1; i>=0; i--) {
 512     if (points[i<<0x1]>-1) {
 513       lrOcc[nCl++] = i;
 514       track.SetClusterIndex(i,points[i<<0x1]);
 515     }
 516   }
 517
 518   if (nCl<3) {
 519     AliError(Form("Cannot estimate momentum of tracks with %d clusters",nCl));
 520     return kFALSE;
 521   }
 522   //
 523   const int lr0   = lrOcc[0];
 524   const int lr1   = lrOcc[nCl/2];
 525   const int lr2   = lrOcc[nCl-1];
 526   //
 527   const SpacePoint& cl0 = fLayer[lr0].Points[ points[lr0<<1] ];
 528   const SpacePoint& cl1 = fLayer[lr1].Points[ points[lr1<<1] ];
 529   const SpacePoint& cl2 = fLayer[lr2].Points[ points[lr2<<1] ];
 530   double cv = Curvature(cl0.XYZ[0],cl0.XYZ[1], cl1.XYZ[0],cl1.XYZ[1], cl2.XYZ[0],cl2.XYZ[1]);
 531
 532   double tgl = (cl2.XYZ[2]-cl0.XYZ[2])/TMath::Sqrt((cl2.XYZ[0]-cl0.XYZ[0])*(cl2.XYZ[0]-cl0.XYZ[0])+(cl2.XYZ[1]-cl0.XYZ[1])*(cl2.XYZ[1]-cl0.XYZ[1]));
 533   //
 534   AliITSUClusterPix* clus = (AliITSUClusterPix*)fClustersTC[ lr0 ]->At( points[lr0<<1] );
 535   AliITSURecoLayer* lr = fITS->GetLayerActive(lr0);
 536   AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
 537   double x = sens->GetXTF() + clus->GetX();
 538   double alp = sens->GetPhiTF();
 539   //  printf("Alp: %f phi: %f\n",alp,phi);
 540   double par[5] = {clus->GetY(),clus->GetZ(),0,tgl,cv};
 541   double cov[15] = {
 542     5*5,
 543     0,  5*5,
 544     0,  0  , 0.7*0.7,
 545     0,  0,   0,       0.7*0.7,
 546     0,  0,   0,       0,      10
 547   };
 548   track.Set(x,alp,par,cov);
 549
 550   return kTRUE;
 551 }
 552
 553 //______________________________________________________________________________
 554 void AliITSUTrackerSA::CandidatesTreeTraversal(vector<Road> &candidates, const int &iD, const int &doubl) {
 555
 556   if ( doubl < 0 ) return;
 557
 558   candidates.back().AddElement(doubl,iD);
 559   const int currentN = candidates.back().N;
 560   for ( size_t iN = 0; iN < fCells[doubl][iD].Neighbours.size(); ++iN ) {
 561     const int currD = doubl - 1 ;
 562     const int neigh = fCells[doubl][iD].Neighbours[iN];
 563
 564     if ( iN > 0 ) {
 565       candidates.push_back(static_cast<Road>(candidates.back()));
 566       candidates.back().N = currentN;
 567     }
 568
 569     CandidatesTreeTraversal(candidates,neigh,currD);
 570   }
 571
 572   fCells[doubl][iD].Level = -1;
 573
 574 }
 575
 576 //______________________________________________________________________________
 577 Double_t AliITSUTrackerSA::RefitTrack(AliExternalTrackParam* trc,
 578               Int_t clInfo[2*AliITSUAux::kMaxLayers],
 579               Double_t rDest, Int_t stopCond)
 580 {
 581   // refit track till radius rDest.
 582   // if stopCond<0 : propagate till last cluster then stop
 583   // if stopCond==0: propagate till last cluster then try to go till limiting rDest, don't mind if fail
 584   // if stopCond>0 : rDest must be reached
 585   //
 586   // The clList should provide the indices of clusters at corresponding layer (as stored in the layer
 587   // TClonesArray, with convention (allowing for up to 2 clusters per layer due to the overlaps):
 588   // if there is a cluster on given layer I, then it should be stored at clInfo[2*I-1]
 589   // if there is an additional cluster on this layer, it goes to clInfo[2*I]
 590   // -1 means no cluster
 591   //
 592   double rCurr = Sqrt(trc->GetX()*trc->GetX() + trc->GetY()*trc->GetY());
 593   int dir,lrStart,lrStop;
 594   //
 595   dir = rCurr<rDest ? 1 : -1;
 596   lrStart = fITS->FindFirstLayerID(rCurr,dir);
 597   lrStop  = fITS->FindLastLayerID(rDest,dir); // lr id before which we have to stop
 598   //
 599   if (lrStop<0 || lrStart<0) AliFatal(Form("Failed to find start(%d) or last(%d) layers. "
 600              "Track from %.3f to %.3f",lrStart,lrStop,rCurr,rDest));
 601   //
 602   int nCl = 0;
 603   for (int i=2*fITS->GetNLayersActive();i--;) {if (clInfo[i]<0) continue; nCl++;}
 604   //
 605   AliExternalTrackParam tmpTr(*trc);
 606   double chi2 = 0;
 607   int iclLr[2],nclLr;
 608   int nclFit = 0;
 609   //
 610   int lrStop1 = lrStop+dir;
 611   for (int ilr=lrStart;ilr!=lrStop1;ilr+=dir) {
 612     AliITSURecoLayer* lr = fITS->GetLayer(ilr);
 613     if ( dir*(rCurr-lr->GetR(dir))>0) continue; // this layer is already passed
 614     int ilrA2,ilrA = lr->GetActiveID();
 615     // passive layer or active w/o hits will be traversed on the way to next cluster
 616     if (!lr->IsActive() || clInfo[ilrA2=(ilrA<<1)]<0) continue;
 617     //
 618     // select the order in which possible 2 clusters (in case of the overlap) will be traversed and fitted
 619     nclLr=0;
 620     if (dir>0) { // clusters are stored in increasing radius order
 621       iclLr[nclLr++]=clInfo[ilrA2++];
 622       if (clInfo[ilrA2]>=0) iclLr[nclLr++]=clInfo[ilrA2];
 623     }
 624     else {
 625       if ( clInfo[ilrA2+1]>=0 ) iclLr[nclLr++]=clInfo[ilrA2+1];
 626       iclLr[nclLr++]=clInfo[ilrA2];
 627     }
 628     //
 629     Bool_t transportedToLayer = kFALSE;
 630     for (int icl=0;icl<nclLr;icl++) {
 631       AliITSUClusterPix* clus =  (AliITSUClusterPix*)lr->GetCluster(iclLr[icl]);
 632       AliITSURecoSens* sens = lr->GetSensorFromID(clus->GetVolumeId());
 633       if (!tmpTr.Rotate(sens->GetPhiTF())) return -1;
 634       //
 635       double xClus = sens->GetXTF()+clus->GetX();
 636       if (!transportedToLayer) {
 637   if (ilr!=lrStart && !TransportToLayerX(&tmpTr,lrStart,ilr,xClus)) return -1; // go to the entrance to the layer
 638   lrStart = ilr;
 639   transportedToLayer = kTRUE;
 640       }
 641       //
 642       if (!PropagateSeed(&tmpTr,xClus,fCurrMass)) return -1;
 643       //
 644       Double_t p[2]={clus->GetY(), clus->GetZ()};
 645       Double_t cov[3]={clus->GetSigmaY2(), clus->GetSigmaYZ(), clus->GetSigmaZ2()};
 646       double chi2cl = tmpTr.GetPredictedChi2(p,cov);
 647       chi2 += chi2cl;
 648       //
 649       if ( !tmpTr.Update(p,cov) ) return -1;
 650       if (++nclFit==nCl && stopCond<0) {
 651   *trc = tmpTr;
 652   return chi2; // it was requested to not propagate after last update
 653       }
 654     }
 655     //
 656   }
 657   // All clusters were succesfully fitted. Even if the track does not reach rDest, this is enough to validate it.
 658   // Still, try to go as close as possible to rDest.
 659   //
 660   if (lrStart!=lrStop) {
 661     if (!TransportToLayer(&tmpTr,lrStart,lrStop)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
 662     if (!GoToExitFromLayer(&tmpTr,fITS->GetLayer(lrStop),dir)) return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
 663   }
 664   // go to the destination radius. Note that here we don't select direction to avoid precision problems
 665   if (!tmpTr.GetXatLabR(rDest,rDest,GetBz(),0) || !PropagateSeed(&tmpTr,rDest,fCurrMass, 100, kFALSE)) {
 666     return (stopCond>0) ? -chi2 : chi2; // rDest was obligatory
 667   }
 668   *trc = tmpTr;
 669
 670   return chi2;
 671 }
 672
 673 //______________________________________________________________________________
 674 Bool_t AliITSUTrackerSA::PropagateSeed(AliExternalTrackParam *seed, Double_t xToGo, Double_t mass, Double_t maxStep, Bool_t matCorr)
 675 {
 676   // propagate seed to given x applying material correction if requested
 677   const Double_t kEpsilon = 1e-5;
 678   Double_t xpos     = seed->GetX();
 679   Int_t dir         = (xpos<xToGo) ? 1:-1;
 680   Double_t xyz0[3],xyz1[3];
 681   //
 682   Bool_t updTime = dir>0 && seed->IsStartedTimeIntegral();
 683   if (matCorr || updTime) seed->GetXYZ(xyz1);   //starting global position
 684   while ( (xToGo-xpos)*dir > kEpsilon){
 685     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 686     Double_t x    = xpos+step;
 687     Double_t bz=GetBz();   // getting the local Bz
 688     if (!seed->PropagateTo(x,bz))  return kFALSE;
 689     double ds = 0;
 690     if (matCorr || updTime) {
 691       xyz0[0]=xyz1[0]; // global pos at the beginning of step
 692       xyz0[1]=xyz1[1];
 693       xyz0[2]=xyz1[2];
 694       seed->GetXYZ(xyz1);    //  // global pos at the end of step
 695       //
 696       if (matCorr) {
 697   Double_t xrho,xx0;
 698   ds = GetMaterialBudget(xyz0,xyz1,xx0,xrho);
 699   if (dir>0) xrho = -xrho; // outward should be negative
 700   if (!seed->CorrectForMeanMaterial(xx0,xrho,mass)) return kFALSE;
 701       }
 702       else { // matCorr is not requested but time integral is
 703   double d0 = xyz1[0]-xyz0[0];
 704   double d1 = xyz1[1]-xyz0[1];
 705   double d2 = xyz1[2]-xyz0[2];
 706   ds = TMath::Sqrt(d0*d0+d1*d1+d2*d2);
 707       }
 708     }
 709     if (updTime) seed->AddTimeStep(ds);
 710     //
 711     xpos = seed->GetX();
 712   }
 713   return kTRUE;
 714 }
 715
 716 //_________________________________________________________________________
 717 Bool_t AliITSUTrackerSA::TransportToLayer(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t rLim)
 718 {
 719   // transport track from layerFrom to the entrance of layerTo or to rLim (if>0), wathever is closer
 720   //
 721   if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
 722   //
 723   int dir = lTo > lFrom ? 1:-1;
 724   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 725   Bool_t checkFirst = kTRUE;
 726   Bool_t limReached = kFALSE;
 727   while(lFrom!=lTo) {
 728     if (lrFr) {
 729       if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
 730       checkFirst = kFALSE;
 731     }
 732     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 733     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 734     //
 735     // go the entrance of the layer, assuming no materials in between
 736     double xToGo = lrTo->GetR(-dir);
 737     if (rLim>0) {
 738       if (dir>0) {
 739   if (rLim<xToGo) {xToGo = rLim; limReached = kTRUE;}
 740       }
 741       else {
 742   if (rLim>xToGo) {xToGo = rLim; limReached = kTRUE;}
 743       }
 744     }
 745     //    double xts = xToGo;
 746     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
 747       //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
 748       //      seed->Print("etp");
 749       return kFALSE;
 750     }
 751     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
 752       //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
 753       //seed->Print("etp");
 754       return kFALSE;
 755     }
 756     lrFr = lrTo;
 757     if (limReached) break;
 758   }
 759   return kTRUE;
 760   //
 761 }
 762
 763 //_________________________________________________________________________
 764 Bool_t AliITSUTrackerSA::TransportToLayerX(AliExternalTrackParam* seed, Int_t lFrom, Int_t lTo, Double_t xStop)
 765 {
 766   // transport track from layerFrom to the entrance of layerTo but do not pass control parameter X
 767   //
 768   if (lTo==lFrom) AliFatal(Form("was called with lFrom=%d lTo=%d",lFrom,lTo));
 769   //
 770   int dir = lTo > lFrom ? 1:-1;
 771   AliITSURecoLayer* lrFr = fITS->GetLayer(lFrom); // this can be 0 when extrapolation from TPC to ITS is requested
 772   Bool_t checkFirst = kTRUE;
 773   while(lFrom!=lTo) {
 774     if (lrFr) {
 775       if (!GoToExitFromLayer(seed,lrFr,dir,checkFirst)) return kFALSE; // go till the end of current layer
 776       checkFirst = kFALSE;
 777     }
 778     AliITSURecoLayer* lrTo =  fITS->GetLayer( (lFrom+=dir) );
 779     if (!lrTo) AliFatal(Form("Layer %d does not exist",lFrom));
 780     //
 781     // go the entrance of the layer, assuming no materials in between
 782     double xToGo = lrTo->GetR(-dir); // R of the entrance to layer
 783     //
 784     //    double xts = xToGo;
 785     if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) {
 786       //      printf("FailHere1: %f %f %d\n",xts,xToGo,dir);
 787       //      seed->Print("etp");
 788       return kFALSE;
 789     }
 790     if ( (dir>0&&xToGo>xStop) || (dir<0&&xToGo<xStop) ) xToGo = xStop;
 791     //
 792 #ifdef _ITSU_DEBUG_
 793     AliDebug(2,Form("go in dir=%d to R=%.4f(X:%.4f)",dir,lrTo->GetR(-dir), xToGo));
 794 #endif
 795     if (!PropagateSeed(seed,xToGo,fCurrMass,100, kFALSE )) {
 796       //printf("FailHere2: %f %f %d\n",xts,xToGo,dir);
 797       //seed->Print("etp");
 798       return kFALSE;
 799     }
 800     lrFr = lrTo;
 801   }
 802   return kTRUE;
 803   //
 804 }
 805
 806 //_________________________________________________________________________
 807 Bool_t AliITSUTrackerSA::GoToExitFromLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 808 {
 809   // go to the exit from lr in direction dir, applying material corrections in steps specific for this layer
 810   // If check is requested, do this only provided the track has not exited the layer already
 811   double xToGo = lr->GetR(dir);
 812   if (check) { // do we need to track till the surface of the current layer ?
 813     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 814     if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // on the surface or outside of the layer
 815     else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // on the surface or outside of the layer
 816   }
 817   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 818   // go via layer to its boundary, applying material correction.
 819   if (!PropagateSeed(seed,xToGo,fCurrMass, lr->GetMaxStep())) return kFALSE;
 820   //
 821   return kTRUE;
 822   //
 823 }
 824
 825 //_________________________________________________________________________
 826 Bool_t AliITSUTrackerSA::GoToEntranceToLayer(AliExternalTrackParam* seed, AliITSURecoLayer* lr, Int_t dir, Bool_t check)
 827 {
 828   // go to the entrance of lr in direction dir, w/o applying material corrections.
 829   // If check is requested, do this only provided the track did not reach the layer already
 830   double xToGo = lr->GetR(-dir);
 831   if (check) { // do we need to track till the surface of the current layer ?
 832     double curR2 = seed->GetX()*seed->GetX() + seed->GetY()*seed->GetY(); // current radius
 833     if      (dir>0) { if (curR2-xToGo*xToGo>-fgkToler) return kTRUE; } // already passed
 834     else if (dir<0) { if (xToGo*xToGo-curR2>-fgkToler) return kTRUE; } // already passed
 835   }
 836   if (!seed->GetXatLabR(xToGo,xToGo,GetBz(),dir)) return kFALSE;
 837   // go via layer to its boundary, applying material correction.
 838   if (!PropagateSeed(seed,xToGo,fCurrMass, 100, kFALSE)) return kFALSE;
 839   return kTRUE;
 840   //
 841 }
 842
 843 //____________________________________________________
 844 Double_t AliITSUTrackerSA::GetMaterialBudget(const double* pnt0,const double* pnt1, double& x2x0, double& rhol) const
 845 {
 846   double par[7];
 847   if (fUseMatLUT && fMatLUT) {
 848     double d = fMatLUT->GetMatBudget(pnt0,pnt1,par);
 849     x2x0 = par[AliITSUMatLUT::kParX2X0];
 850     rhol = par[AliITSUMatLUT::kParRhoL];
 851     return d;
 852   }
 853   else {
 854     MeanMaterialBudget(pnt0,pnt1,par);
 855     x2x0 = par[1];
 856     rhol = par[0]*par[4];
 857     return par[4];
 858   }
 859 }
 860
 861 //____________________________________________________________________
 862 Double_t AliITSUTrackerSA::Curvature(Double_t x1,Double_t y1,Double_t x2,Double_t y2,Double_t x3,Double_t y3) {
 863
 864   //calculates the curvature of track
 865   Double_t den = (x3-x1)*(y2-y1)-(x2-x1)*(y3-y1);
 866   if(den*den<1e-32) return 0.;
 867   Double_t a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 868   if((y2-y1)*(y2-y1)<1e-32) return 0.;
 869   Double_t b = -(x2*x2-x1*x1+y2*y2-y1*y1+a*(x2-x1))/(y2-y1);
 870   Double_t c = -x1*x1-y1*y1-a*x1-b*y1;
 871   Double_t xc= -a/2.;
 872
 873   if((a*a+b*b-4*c)<0) return 0.;
 874   Double_t rad = TMath::Sqrt(a*a+b*b-4*c)/2.;
 875   if(rad*rad<1e-32) return 1e16;
 876
 877   if((x1>0 && y1>0 && x1<xc)) rad*=-1;
 878   if((x1<0 && y1>0 && x1<xc)) rad*=-1;
 879     //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
 880     // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
 881
 882   return 1/rad;
 883
 884 }
 885