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