ITS/UPGRADE/AliITSUCATracker.cxx

   1 //-*- Mode: C++ -*-
   2 // **************************************************************************
   3 // This file is property of and copyright by the ALICE ITSU Project         *
   4 // ALICE Experiment at CERN, All rights reserved.                           *
   5 //                                                                          *
   6 // Primary Author: Maximiliano Puccio <maximiliano.puccio@cern.ch>          *
   7 //                 for the ITS Upgrade project                              *
   8 //                                                                          *
   9 // Permission to use, copy, modify and distribute this software and its     *
  10 // documentation strictly for non-commercial purposes is hereby granted     *
  11 // without fee, provided that the above copyright notice appears in all     *
  12 // copies and that both the copyright notice and this permission notice     *
  13 // appear in the supporting documentation. The authors make no claims       *
  14 // about the suitability of this software for any purpose. It is            *
  15 // provided "as is" without express or implied warranty.                    *
  16 //                                                                          *
  17 //***************************************************************************
  18
  19 #include "AliITSUCATracker.h"
  20
  21 // STD
  22 #include <algorithm>
  23 #include <cassert>
  24 // ROOT
  25 #include <TBranch.h>
  26 #include <TMath.h>
  27 #include <TTree.h>
  28 #include <Riostream.h>
  29 // ALIROOT
  30 #include "AliESDEvent.h"
  31 #include "AliESDtrack.h"
  32 #include "AliLog.h"
  33 // ALIROOT ITSU
  34 #include "AliITSUAux.h"
  35 #include "AliITSUClusterPix.h"
  36 #include "AliITSURecoDet.h"
  37 #include "AliITSUReconstructor.h"
  38 #include "AliITSUTrackCooked.h"
  39
  40 using TMath::Abs;
  41 using TMath::Sort;
  42 using TMath::Sqrt;
  43 using std::sort;
  44 using std::cout;
  45 using std::endl;
  46 using std::flush;
  47
  48 ClassImp(AliITSUCATracker)
  49
  50
  51 // tolerance for layer on-surface check
  52 const Double_t AliITSUCATracker::fgkChi2Cut =  600.f;
  53 const int AliITSUCATracker::fgkNumberOfIterations =  2;
  54 const float AliITSUCATracker::fgkR[7] = {2.33959,3.14076,3.91924,19.6213,24.5597,34.388,39.3329};
  55 //
  56 const float kmaxDCAxy[5] = {0.05f,0.04f,0.05f,0.2f,0.4f};
  57 const float kmaxDCAz[5] = {0.2f,0.4f,0.5f,0.6f,3.f};
  58 const float kmaxDN[4] = {0.002f,0.009f,0.002f,0.005f};
  59 const float kmaxDP[4] = {0.008f,0.0025f,0.003f,0.0035f};
  60 const float kmaxDZ[6] = {0.1f,0.1f,0.3f,0.3f,0.3f,0.3f};
  61 const float kDoublTanL = 0.025;
  62 const float kDoublPhi = 0.14;
  63
  64 const float kmaxDCAxy1[5] = /*{1.f,0.5,0.5,1.7,3.};/*/{1.f,0.4f,0.4f,1.5f,3.f};
  65 const float kmaxDCAz1[5] = /*{2.f,0.8,0.8,3.,5.};/*/{1.f,0.4f,0.4f,1.5f,3.f};
  66 const float kmaxDN1[4] = /*{0.006f,0.0045f,0.01f,0.04f};/*/{0.005f,0.0035f,0.009f,0.03f};
  67 const float kmaxDP1[4] = /*{0.04f,0.01f,0.012f,0.014f};/*/{0.02f,0.005f,0.006f,0.007f};
  68 const float kmaxDZ1[6] = /*{1.5f,1.5f,2.f,2.f,2.f,2.f};/*/{1.f,1.f,1.5f,1.5f,1.5f,1.5f};
  69 const float kDoublTanL1 = /*0.12f;/*/0.05f;
  70 const float kDoublPhi1 = /*0.4f;/*/0.2f;
  71
  72 //__________________________________________________________________________________________________
  73 static inline float invsqrt(float _x)
  74 {
  75   //
  76   // The function calculates fast inverse sqrt. Google for 0x5f3759df.
  77   // Credits to ALICE HLT Project
  78   //
  79
  80   union { float f; int i; } x = { _x };
  81   const float xhalf = 0.5f * x.f;
  82   x.i = 0x5f3759df - ( x.i >> 1 );
  83   x.f = x.f * ( 1.5f - xhalf * x.f * x.f );
  84   return x.f;
  85 }
  86
  87 //__________________________________________________________________________________________________
  88 static inline float Curvature(float x1, float y1, float x2, float y2, float x3, float y3)
  89 {
  90   //
  91   // Initial approximation of the track curvature
  92   //
  93 //  return - 2.f * ((x2 - x1) * (y3 - y2) - (x3 - x2) * (y2 - y1))
  94 //               * invsqrt(((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1)) *
  95 //                         ((x2 - x3) * (x2 - x3) + (y2 - y3) * (y2 - y3)) *
  96 //                         ((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1)));
  97
  98   //calculates the curvature of track
  99   float den = (x3 - x1) * (y2 - y1) - (x2 - x1) * (y3 - y1);
 100   if(den * den < 1e-32) return 0.f;
 101   float a = ((y3-y1)*(x2*x2+y2*y2-x1*x1-y1*y1)-(y2-y1)*(x3*x3+y3*y3-x1*x1-y1*y1))/den;
 102   if((y2 - y1) * (y2 - y1) < 1e-32) return 0.f;
 103   float b = -(x2 * x2 - x1 * x1 + y2 * y2 - y1 * y1 + a * (x2 - x1)) / (y2 - y1);
 104   float c = -x1 * x1 - y1 * y1 - a * x1 - b * y1;
 105   float xc = -a / 2.f;
 106
 107   if((a * a + b * b - 4 * c) < 0) return 0.f;
 108   float rad = TMath::Sqrt(a * a + b * b - 4 * c) / 2.f;
 109   if(rad * rad < 1e-32) return 1e16;
 110
 111   if((x1 > 0.f && y1 > 0.f && x1 < xc)) rad *= -1.f;
 112   if((x1 < 0.f && y1 > 0.f && x1 < xc)) rad *= -1.f;
 113   //  if((x1<0 && y1<0 && x1<xc)) rad*=-1;
 114   // if((x1>0 && y1<0 && x1<xc)) rad*=-1;
 115
 116   return 1.f/rad;
 117
 118 }
 119
 120 //__________________________________________________________________________________________________
 121 static inline float TanLambda(float x1, float y1, float x2, float y2, float z1, float z2)
 122 {
 123   //
 124   // Initial approximation of the tangent of the track dip angle
 125   //
 126   return (z1 - z2) * invsqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
 127 }
 128
 129 //__________________________________________________________________________________________________
 130 //static inline float XCenterOfCurvature(float x1, float y1, float x2, float y2, float x3, float y3)
 131 //{
 132 //    //
 133 //    // Initial approximation of the x-coordinate of the center of curvature
 134 //    //
 135 //
 136 //  const float k1 = (y2 - y1) / (x2 - x1), k2 = (y3 - y2) / (x3 - x2);
 137 //  return TMath::Abs(k2 - k1) > kAlmost0 ?
 138 //    0.5f * (k1 * k2 * (y1 - y3) + k2 * (x1 + x2) - k1 * (x2 + x3)) / (k2 - k1) : 1e12f;
 139 //}
 140
 141 //__________________________________________________________________________________________________
 142 static inline bool CompareAngles(float alpha, float beta, float tolerance)
 143 {
 144         const float delta = TMath::Abs(alpha - beta);
 145         return (delta < tolerance || TMath::Abs(delta - TMath::TwoPi()) < tolerance);
 146 }
 147
 148 //__________________________________________________________________________________________________
 149 AliITSUCATracker::AliITSUCATracker(AliITSUReconstructor* rec) : AliITSUTrackerGlo(rec)
 150 #ifdef _TUNING_
 151 ,fGood(0)
 152 ,fTan(NULL)
 153 ,fTanF(NULL)
 154 ,fPhiF(NULL)
 155 ,fNEntries(NULL)
 156 #endif
 157 ,fLayer()
 158 ,fUsedClusters()
 159 ,fChi2Cut(fgkChi2Cut)
 160 ,fPhiCut(1)
 161 ,fZCut(0.5f)
 162 ,fCandidates()
 163 ,fSAonly(kTRUE)
 164 ,fCPhi()
 165 ,fCDTanL()
 166 ,fCDPhi()
 167 ,fCZ()
 168 ,fCDCAz()
 169 ,fCDCAxy()
 170 ,fCDN()
 171 ,fCDP()
 172 ,fCDZ()
 173 {
 174   // This default constructor needs to be provided
 175   for (int i = 0; i < 4; ++i)
 176   {
 177     fCandidates[i] = new TClonesArray("AliITSUTrackCooked",100000);
 178   }
 179
 180 #ifdef _TUNING_
 181   for (int i = 0; i < 6; ++i)
 182   {
 183     fGDZ[i] = new TH1F(Form("DGZ%i",i),Form("DZ%i;#Deltaz;Entries",i),500,0.f,5.f);
 184     fGDXY[i] = new TH1F(Form("DGPhi%i",i),Form("#Delta#Phi%i;#DeltaPhi;Entries",i),500,0.f,TMath::TwoPi());
 185     fFDZ[i] = new TH1F(Form("DFZ%i",i),Form("DZ%i;#Deltaz;Entries",i),500,0.f,5.f);
 186     fFDXY[i] = new TH1F(Form("DFPhi%i",i),Form("#Delta#Phi%i;#Delta#Phi;Entries",i),500,0.f,TMath::TwoPi());
 187   }
 188   for (int i = 0; i < 5; ++i)
 189   {
 190     fGDCAZ[i] = new TH1F(Form("DCAGZ%i",i),Form("DCAZ%i;#Deltaz;Entries",i),500,0.f,5.f);
 191     fGDCAXY[i] = new TH1F(Form("DCAGXY%i",i),Form("DCAXY%i;#Deltar;Entries",i),500,0.f,5.f);
 192     fFDCAZ[i] = new TH1F(Form("DCAFZ%i",i),Form("DCAZ%i;#Deltaz;Entries",i),500,0.f,5.f);
 193     fFDCAXY[i] = new TH1F(Form("DCAFXY%i",i),Form("DCAXY%i;#Deltar;Entries",i),500,0.f,5.f);
 194   }
 195   for(int i = 0; i < 4; ++i)
 196   {
 197     fGoodCombChi2[i] = new TH1F(Form("goodcombchi2%i",i),Form("%i;#chi^{2};Entries",i),2000,0,0.1);
 198     fFakeCombChi2[i] = new TH1F(Form("fakecombchi2%i",i),Form("%i;#chi^{2};Entries",i),2000,0,0.1);
 199     fGoodCombN[i] = new TH1F(Form("goodcombn%i",i),Form("%i;#Deltan;Entries",i),300,0.f,0.03f);
 200     fFakeCombN[i] = new TH1F(Form("fakecombn%i",i),Form("%i;#Deltan;Entries",i),300,0.f,0.03f);
 201   }
 202   fGoodCombChi2[4] = new TH1F("goodcomb4",";#chi^{2};Entries",200,0,500);
 203   fFakeCombChi2[4] = new TH1F("fakecomb4",";#chi^{2};Entries",200,0,500);
 204   fTan = new TH1F("tan","tan",2500,0,0.5);
 205   fTanF = new TH1F("tanF","tanF",2500,0,0.5);
 206   fPhi = new TH1F("phi","phi",2500,0,TMath::Pi());
 207   fPhiF = new TH1F("phi","phiF",2500,0,TMath::Pi());
 208   fNEntries = new TH1F("nentries","nentries",2001,-0.5,2000.5);
 209 #endif
 210 }
 211
 212 //__________________________________________________________________________________________________
 213 AliITSUCATracker::~AliITSUCATracker()
 214 {
 215   // d-tor
 216    for (int i = 0; i < 4; ++i)
 217   {
 218     if (fCandidates[i])
 219       delete fCandidates[i];
 220   }
 221 //#ifdef _TUNING_
 222 //  // Just cut and paste from the constructor
 223 //  for (int i = 0; i < 6; ++i)
 224 //  {
 225 //    delete fGDZ[i];
 226 //    delete fGDXY[i];
 227 //    delete fFDZ[i];
 228 //    delete fFDXY[i];
 229 //  }
 230 //  for (int i = 0; i < 5; ++i)
 231 //  {
 232 //    delete fGDCAZ[i];
 233 //    delete fGDCAXY[i];
 234 //    delete fFDCAZ[i];
 235 //    delete fFDCAXY[i];
 236 //  }
 237 //  for(int i = 0; i < 4; ++i)
 238 //  {
 239 //    delete fGoodCombChi2[i];
 240 //    delete fFakeCombChi2[i];
 241 //    delete fGoodCombN[i];
 242 //    delete fFakeCombN[i];
 243 //  }
 244 //  delete fGoodCombChi2[4];
 245 //  delete fFakeCombChi2[4];
 246 //  delete fTan;
 247 //  delete fTanF;
 248 //  delete fPhi;
 249 //  delete fPhiF;
 250 //  delete fNEntries;
 251 //#endif
 252 }
 253
 254 //__________________________________________________________________________________________________
 255 bool AliITSUCATracker::CellParams(int l, ClsInfo_t* c1, ClsInfo_t* c2, ClsInfo_t* c3,
 256                                   float &curv, float n[3])
 257 {
 258   // Calculation of cell params and filtering using a DCA cut wrt beam line position.
 259   // The hit are mapped on a paraboloid space: there a circle is described as plane.
 260   // The parameter n of the cells is the normal vector to the plane describing the circle in the
 261   // paraboloid.
 262
 263   // Mapping the hits
 264   const float mHit0[3] = {c1->x, c1->y, c1->r * c1->r};
 265   const float mHit1[3] = {c2->x, c2->y, c2->r * c2->r};
 266   const float mHit2[3] = {c3->x, c3->y, c3->r * c3->r};
 267   // Computing the deltas
 268   const float mD10[3] = {mHit1[0] - mHit0[0],mHit1[1] - mHit0[1],mHit1[2] - mHit0[2]};
 269   const float mD20[3] = {mHit2[0] - mHit0[0],mHit2[1] - mHit0[1],mHit2[2] - mHit0[2]};
 270   // External product of the deltas -> n
 271   n[0] = (mD10[1] * mD20[2]) - (mD10[2] * mD20[1]);
 272   n[1] = (mD10[2] * mD20[0]) - (mD10[0] * mD20[2]);
 273   n[2] = (mD10[0] * mD20[1]) - (mD10[1] * mD20[0]);
 274   // Normalisation
 275   const float norm = TMath::Sqrt((n[0] * n[0]) + (n[1] * n[1]) + (n[2] * n[2]));
 276   if (norm < 1e-20f || fabs(n[2]) < 1e-20f)
 277     return false;
 278   n[0] /= norm;
 279   n[1] /= norm;
 280   n[2] /= norm;
 281   // Center of the circle
 282   const float c[2] = {-0.5f * n[0] / n[2], -0.5f * n[1] / n[2]};
 283   // Constant
 284   const float k = - n[0] * mHit1[0] - n[1] * mHit1[1] - n[2] * mHit1[2];
 285   // Radius of the circle
 286   curv = TMath::Sqrt((1.f - n[2] * n[2] - 4.f * k * n[2]) / (4.f * n[2] * n[2]));
 287   // Distance of closest approach to the beam line
 288   const float dca = fabs(curv - sqrt(c[0] * c[0] + c[1] * c[1]));
 289   // Cut on the DCA
 290   if (dca > fCDCAxy[l]) {
 291     return false;
 292   }
 293 #ifdef _TUNING_
 294   if (fGood) {
 295     fGDCAXY[l]->Fill(dca);
 296   } else {
 297     fFDCAXY[l]->Fill(dca);
 298   }
 299 #endif
 300
 301   curv = 1.f / curv;
 302   return true;
 303 }
 304
 305 //__________________________________________________________________________________________________
 306 void AliITSUCATracker::CellsTreeTraversal(vector<AliITSUCARoad> &roads,
 307                                           const int &iD, const int &doubl)
 308 {
 309
 310   // Each cells contains a list of neighbours. Each neighbour has presumably other neighbours.
 311   // This chain of neighbours is, as a matter of fact, a tree and this function goes recursively
 312   // through it. This function implements a depth first tree browsing.
 313
 314   // End of the road
 315   if (doubl < 0) return;
 316
 317   // [1] add current cell to current cell
 318   roads.back().AddElement(doubl,iD);
 319   // We want the right number of elements in the roads also in the case of multiple neighbours
 320   const int currentN = roads.back().N;
 321
 322   // [2] loop on the neighbours of the current cell
 323   for (size_t iN = 0; iN < fCells[doubl][iD].NumberOfNeighbours(); ++iN)
 324   {
 325     const int currD = doubl - 1;
 326     const int neigh = fCells[doubl][iD](iN);
 327
 328     // [3] for each neighbour one road
 329     if (iN > 0)
 330     {
 331       roads.push_back(roads.back());
 332       roads.back().N = currentN;
 333     }
 334     // [4] play this game again until the end of the road
 335     CellsTreeTraversal(roads,neigh,currD);
 336   }
 337
 338   fCells[doubl][iD].SetLevel(0u); // Level = -1
 339 }
 340
 341 //__________________________________________________________________________________________________
 342 Int_t AliITSUCATracker::Clusters2Tracks(AliESDEvent *event)
 343 {
 344   // This is the main tracking function
 345   // The clusters must already be loaded
 346
 347   if (!fSAonly) {
 348     AliITSUTrackerGlo::Clusters2Tracks(event);
 349     for (int iL = 0; iL < 7; ++iL) {
 350       for (int iC = 0; iC < fLayer[iL].GetNClusters(); ++iC)
 351         if (fLayer[iL].GetClusterUnSorted(iC)->IsClusterUsed())
 352           fUsedClusters[iL][iC] = true;
 353     }
 354   }
 355
 356   int ntrk = 0, ngood = 0;
 357   for (int iteration = 0; iteration < fgkNumberOfIterations; ++iteration)
 358   {
 359
 360     fCandidates[0]->Clear();
 361     fCandidates[1]->Clear();
 362     fCandidates[2]->Clear();
 363     fCandidates[3]->Clear();
 364
 365     MakeCells(iteration);
 366     FindTracksCA(iteration);
 367
 368     for (int iL = 3; iL >= 0; --iL)
 369     {
 370       const int nCand = fCandidates[iL]->GetEntries();
 371       int index[nCand];
 372       float chi2[nCand];
 373
 374       for (int iC = 0; iC < nCand; ++iC)
 375       {
 376         AliITSUTrackCooked *tr = (AliITSUTrackCooked*)fCandidates[iL]->At(iC);
 377         chi2[iC] = tr->GetChi2();//(RefitTrack(tr,clInfo,0.,-1));
 378         index[iC] = iC;
 379         CookLabel(tr,0.f);
 380       }
 381
 382       TMath::Sort(nCand,chi2,index,false);
 383
 384       for (int iUC = 0; iUC < nCand; ++iUC)
 385       {
 386         const int iC = index[iUC];
 387         if (chi2[iC] < 0.f)
 388         {
 389           continue;
 390         }
 391
 392         AliITSUTrackCooked *tr = (AliITSUTrackCooked*)fCandidates[iL]->At(iC);
 393
 394         bool sharingCluster = false;
 395         for (int k = 0; k < tr->GetNumberOfClusters(); ++k)
 396         {
 397           const int layer = (tr->GetClusterIndex(k) & 0xf0000000) >> 28;
 398           const int idx = (tr->GetClusterIndex(k) & 0x0fffffff);
 399           if (fUsedClusters[layer][idx])
 400           {
 401             sharingCluster = true;
 402             break;
 403           }
 404         }
 405
 406         if (sharingCluster)
 407           continue;
 408
 409         for (int k = 0; k < tr->GetNumberOfClusters(); ++k)
 410         {
 411           const int layer = (tr->GetClusterIndex(k) & 0xf0000000) >> 28;
 412           const int idx = (tr->GetClusterIndex(k) & 0x0fffffff);
 413           fUsedClusters[layer][idx] = true;
 414         }
 415
 416         AliESDtrack outTrack;
 417         CookLabel(tr,0.f);
 418         ntrk++;
 419         if(tr->GetLabel() >= 0)
 420         {
 421           ngood++;
 422 #ifdef _TUNING_
 423           fGoodCombChi2[4]->Fill(chi2[iC] / (4 + iL));
 424         }
 425         else
 426         {
 427           fFakeCombChi2[4]->Fill(chi2[iC] / (4 + iL));
 428 #endif
 429         }
 430
 431         outTrack.UpdateTrackParams(tr,AliESDtrack::kITSin);
 432         outTrack.SetLabel(tr->GetLabel());
 433         event->AddTrack(&outTrack);
 434       }
 435     }
 436   }
 437   Info("Clusters2Tracks","Reconstructed tracks: %d",ntrk);
 438   if (ntrk)
 439     Info("Clusters2Tracks","Good tracks/reconstructed: %f",Float_t(ngood)/ntrk);
 440   //
 441   return 0;
 442 }
 443
 444 //__________________________________________________________________________________________________
 445 void AliITSUCATracker::FindTracksCA(int iteration)
 446 {
 447   // Main pattern recognition routine. It has 4 steps (planning to split in different methods)
 448   // 1. Tracklet finding (using vertex position)
 449   // 2. Tracklet association, cells building
 450   // 3. Handshake between neighbour cells
 451   // 4. Candidates ("roads") finding and fitting
 452
 453   // Road finding and fitting. The routine starts from cells with level 5 since they are the head
 454   // of a full road (candidate with 7 points). Minimum level is 2, so candidates at the end have
 455   // at least 4 points.
 456   const int itLevelLimit[3] = {4, 4, 1};
 457   for (int level = 5; level > itLevelLimit[iteration]; --level) {
 458     vector<AliITSUCARoad> roads;
 459     // Road finding. For each cell at level $(level) a loop on their neighbours to start building
 460     // the roads.
 461     for (int iCL = 4; iCL >= level - 1; --iCL) {
 462       for (size_t iCell = 0; iCell < fCells[iCL].size(); ++iCell) {
 463         if (fCells[iCL][iCell].GetLevel() != level)
 464         {
 465           continue;
 466         }
 467         // [1] Add current cell to road
 468         roads.push_back(AliITSUCARoad(iCL,iCell));
 469         // [2] Loop on current cell neighbours
 470         for(size_t iN = 0; iN < fCells[iCL][iCell].NumberOfNeighbours(); ++iN) {
 471           const int currD = iCL - 1;
 472           const int neigh = fCells[iCL][iCell](iN);
 473           // [3] if more than one neighbour => more than one road, one road for each neighbour
 474           if(iN > 0)
 475           {
 476             roads.push_back(AliITSUCARoad(iCL,iCell));
 477           }
 478           // [4] Essentially the neighbour became the current cell and then go to [1]
 479           CellsTreeTraversal(roads,neigh,currD);
 480         }
 481         fCells[iCL][iCell].SetLevel(0u); // Level = -1
 482       }
 483     }
 484
 485     // Roads fitting
 486     for (size_t iR = 0; iR < roads.size(); ++iR)
 487     {
 488       if (roads[iR].N != level)
 489         continue;
 490       AliITSUTrackCooked tr;
 491       int first = -1,last = -1;
 492       ClsInfo_t *cl0 = 0x0,*cl1 = 0x0,*cl2 = 0x0;
 493       for(int i = 0; i < 5; ++i)
 494       {
 495         if (roads[iR][i] < 0)
 496           continue;
 497
 498         if (first < 0)
 499         {
 500           cl0 = fLayer[i][fCells[i][roads[iR][i]].x()];
 501           tr.SetClusterIndex(i,fLayer[i][fCells[i][roads[iR][i]].x()]->index);
 502           tr.SetClusterIndex(i + 1,fLayer[i + 1][fCells[i][roads[iR][i]].y()]->index);
 503           first = i;
 504         }
 505         tr.SetClusterIndex(i + 2,fLayer[i + 2][fCells[i][roads[iR][i]].z()]->index);
 506         last = i;
 507       }
 508       AliITSUClusterPix* c = fLayer[last + 2].GetClusterSorted(fCells[last][roads[iR][last]].z());
 509       cl2 = fLayer[last + 2][fCells[last][roads[iR][last]].z()];
 510       first = (last + first) / 2;
 511       cl1 = fLayer[first + 1][fCells[first][roads[iR][first]].y()];
 512       // Init track parameters
 513       double cv = Curvature(cl0->x,cl0->y,cl1->x,cl1->y,cl2->x,cl2->y);
 514       double tgl = TanLambda(cl0->x,cl0->y,cl2->x,cl2->y,cl0->z,cl2->z);
 515
 516       AliITSUCATrackingStation::ITSDetInfo_t det = fLayer[last + 2].GetDetInfo(cl2->detid);
 517       double x = det.xTF + c->GetX(); // I'd like to avoit using AliITSUClusterPix...
 518       double alp = det.phiTF;
 519       double par[5] = {c->GetY(),c->GetZ(),0,tgl,cv};
 520       double cov[15] = {
 521         5*5,
 522         0,  5*5,
 523         0,  0  , 0.7*0.7,
 524         0,  0,   0,       0.7*0.7,
 525         0,  0,   0,       0,       10
 526       };
 527       tr.Set(x,alp,par,cov);
 528       AliITSUTrackCooked tt(tr);
 529       tt.ResetClusters();
 530       if (RefitAt(2.1, &tt, &tr))
 531         new((*fCandidates[level - 2])[fCandidates[level - 2]->GetEntriesFast()]) AliITSUTrackCooked(tt);
 532     }
 533   }
 534 }
 535
 536 //__________________________________________________________________________________________________
 537 inline AliCluster* AliITSUCATracker::GetCluster(Int_t index) const
 538 {
 539         const Int_t l=(index & 0xf0000000) >> 28;
 540         const Int_t c=(index & 0x0fffffff);
 541         return (AliCluster*)fLayer[l].GetClusterUnSorted(c) ;
 542 }
 543
 544 //__________________________________________________________________________________________________
 545 Double_t AliITSUCATracker::GetMaterialBudget(const double* pnt0,const double* pnt1, double& x2x0,
 546                                                    double& rhol) const
 547 {
 548   double par[7];
 549   if (fUseMatLUT && fMatLUT) {
 550     double d = fMatLUT->GetMatBudget(pnt0,pnt1,par);
 551     x2x0 = par[AliITSUMatLUT::kParX2X0];
 552     rhol = par[AliITSUMatLUT::kParRhoL];
 553     return d;
 554   }
 555   else {
 556     MeanMaterialBudget(pnt0,pnt1,par);
 557     x2x0 = par[1];
 558     rhol = par[0]*par[4];
 559     return par[4];
 560   }
 561 }
 562
 563 //__________________________________________________________________________________________________
 564 Int_t AliITSUCATracker::LoadClusters(TTree *cluTree)
 565 {
 566   // This function reads the ITSU clusters from the tree,
 567   // sort them, distribute over the internal tracker arrays, etc
 568
 569   AliITSUTrackerGlo::LoadClusters(cluTree); // === fITS->LoadClusters(cluTree);
 570   if (fSAonly) fITS->ProcessClusters();
 571
 572   // I consider a single vertex event for the moment.
 573   //TODO: pile-up (trivial here), fast reco of primary vertices (not so trivial)
 574   double vertex[3] = {GetX(),GetY(),GetZ()};
 575   for (int iL = 0; iL < 7; ++iL) {
 576     fLayer[iL].Init(fITS->GetLayerActive(iL),fITS->GetGeom());
 577     AliVertex v(vertex,1,1);
 578     fLayer[iL].SortClusters(&v);
 579     fUsedClusters[iL].resize(fLayer[iL].GetNClusters(),false);
 580   }
 581   return 0;
 582 }
 583
 584 //__________________________________________________________________________________________________
 585 void AliITSUCATracker::MakeCells(int iteration)
 586 {
 587 #ifdef _TUNING_
 588   unsigned int numberOfGoodDoublets = 0, totalNumberOfDoublets = 0;
 589   unsigned int numberOfGoodCells = 0, totalNumberOfCells = 0;
 590   unsigned int cellsCombiningSuccesses = 0, cellsWrongCombinations = 0;
 591 #endif
 592
 593   SetCuts(iteration);
 594   if (iteration >= 1) {
 595 #ifdef _TUNING_
 596     ResetHistos();
 597 #endif
 598     for (int i = 0; i < 5; ++i)
 599       fCells[i].clear();
 600     for (int i = 0; i < 6; ++i)
 601       fDoublets[i].clear();
 602   }
 603
 604   // Trick to speed up the navigation of the doublets array. The lookup table is build like:
 605   // dLUT[l][i] = n;
 606   // where n is the index inside fDoublets[l+1] of the first doublets that uses the point
 607   // fLayer[l+1][i]
 608   vector<int> dLUT[5];
 609   for (int iL = 0; iL < 6; ++iL) {
 610     if (iL < 5) {
 611       dLUT[iL].resize(fLayer[iL + 1].GetNClusters(),-1);
 612     }
 613     for (int iC = 0; iC < fLayer[iL].GetNClusters(); ++iC) {
 614       ClsInfo_t* cls = fLayer[iL].GetClusterInfo(iC);
 615       if (fUsedClusters[iL][cls->index]) {
 616         continue;
 617       }
 618       const float tanL = (cls->z - GetZ()) / cls->r;
 619       const float extz = tanL * (fgkR[iL + 1] - cls->r) + cls->z;
 620       const int nClust = fLayer[iL + 1].SelectClusters(extz - 2 * fCZ, extz + 2 * fCZ,
 621                                                        cls->phi - fCPhi, cls->phi + fCPhi);
 622       bool first = true;
 623
 624       for (int iC2 = 0; iC2 < nClust; ++iC2) {
 625         const int iD2 = fLayer[iL + 1].GetNextClusterInfoID();
 626         ClsInfo_t* cls2 = fLayer[iL + 1].GetClusterInfo(iD2);
 627         if (fUsedClusters[iL + 1][cls2->index]) {
 628           continue;
 629         }
 630         const float dz = tanL * (cls2->r - cls->r) + cls->z - cls2->z;
 631         if (TMath::Abs(dz) < fCDZ[iL] && CompareAngles(cls->phi, cls2->phi, fCPhi)) {
 632           if (first && iL > 0) {
 633             dLUT[iL - 1][iC] = fDoublets[iL].size();
 634             first = false;
 635           }
 636           const float dTanL = (cls->z - cls2->z) / (cls->r - cls2->r);
 637           const float phi = TMath::ATan2(cls->y - cls2->y, cls->x - cls2->x);
 638           fDoublets[iL].push_back(Doublets(iC,iD2,dTanL,phi));
 639 #ifdef _TUNING_
 640           if (fLayer[iL].GetClusterSorted(iC)->GetLabel(0) ==
 641               fLayer[iL + 1].GetClusterSorted(iD2)->GetLabel(0) &&
 642               fLayer[iL].GetClusterSorted(iC)->GetLabel(0) > 0) {
 643             numberOfGoodDoublets++;
 644             fGDZ[iL]->Fill(dz);
 645             fGDXY[iL]->Fill(fabs(cls->phi - cls2->phi));
 646           } else {
 647             fFDZ[iL]->Fill(dz);
 648             fFDXY[iL]->Fill(fabs(cls->phi - cls2->phi));
 649           }
 650           totalNumberOfDoublets++;
 651 #endif
 652         }
 653       }
 654       fLayer[iL + 1].ResetFoundIterator();
 655     }
 656   }
 657
 658   // Trick to speed up the navigation of the cells array. The lookup table is build like:
 659   // tLUT[l][i] = n;
 660   // where n is the index inside fCells[l+1] of the first cells that uses the doublet
 661   // fDoublets[l+1][i]
 662   vector<int> tLUT[4];
 663   for (int iD = 0; iD < 5; ++iD)
 664   {
 665     if (iD < 4) {
 666       tLUT[iD].resize(fDoublets[iD + 1].size(),0);
 667     }
 668     for (size_t iD0 = 0; iD0 < fDoublets[iD].size(); ++iD0)
 669     {
 670       const int idx = fDoublets[iD][iD0].y;
 671       bool first = true;
 672       if (dLUT[iD][idx] == -1) continue;
 673       for (size_t iD1 = dLUT[iD][idx]; idx == fDoublets[iD + 1][iD1].x;++iD1)
 674       {
 675         if (TMath::Abs(fDoublets[iD][iD0].tanL - fDoublets[iD + 1][iD1].tanL) < fCDTanL &&
 676             TMath::Abs(fDoublets[iD][iD0].phi - fDoublets[iD + 1][iD1].phi) < fCDPhi) {
 677           const float tan = 0.5f * (fDoublets[iD][iD0].tanL + fDoublets[iD + 1][iD1].tanL);
 678           const float extz = -tan * fLayer[iD][fDoublets[iD][iD0].x]->r +
 679                               fLayer[iD][fDoublets[iD][iD0].x]->z;
 680           if (fabs(extz - GetZ()) < fCDCAz[iD]) {
 681 #ifdef _TUNING_
 682             fGood = (fLayer[iD].GetClusterSorted(fDoublets[iD][iD0].x)->GetLabel(0) ==
 683                      fLayer[iD + 1].GetClusterSorted(fDoublets[iD][iD0].y)->GetLabel(0) &&
 684                      fLayer[iD].GetClusterSorted(fDoublets[iD][iD0].x)->GetLabel(0) ==
 685                      fLayer[iD + 2].GetClusterSorted(fDoublets[iD + 1][iD1].y)->GetLabel(0) &&
 686                      fLayer[iD].GetClusterSorted(fDoublets[iD][iD0].x)->GetLabel(0) > 0);
 687 #endif
 688             float curv, n[3];
 689             if (CellParams(iD,fLayer[iD][fDoublets[iD][iD0].x],fLayer[iD + 1][fDoublets[iD][iD0].y],
 690                            fLayer[iD + 2][fDoublets[iD + 1][iD1].y],curv,n)) {
 691               if (first && iD > 0) {
 692                 tLUT[iD - 1][iD0] = fCells[iD].size();
 693                 first = false;
 694               }
 695               fCells[iD].push_back(AliITSUCACell(fDoublets[iD][iD0].x,fDoublets[iD][iD0].y,
 696                                                  fDoublets[iD + 1][iD1].y,iD0,iD1,curv,n));
 697 #ifdef _TUNING_
 698               if (fGood) {
 699                 fTan->Fill(TMath::Abs(fDoublets[iD][iD0].tanL - fDoublets[iD + 1][iD1].tanL));
 700                 fPhi->Fill(TMath::Abs(fDoublets[iD][iD0].phi - fDoublets[iD + 1][iD1].phi));
 701                 fGDCAZ[iD]->Fill(fabs(extz-GetZ()));
 702                 numberOfGoodCells++;
 703               } else {
 704                 fTanF->Fill(TMath::Abs(fDoublets[iD][iD0].tanL - fDoublets[iD + 1][iD1].tanL));
 705                 fPhiF->Fill(TMath::Abs(fDoublets[iD][iD0].phi - fDoublets[iD + 1][iD1].phi));
 706                 fFDCAZ[iD]->Fill(fabs(extz - GetZ()));
 707               }
 708               totalNumberOfCells++;
 709 #endif
 710             }
 711           }
 712         }
 713       }
 714     }
 715   }
 716
 717   // Adjacent cells: cells that share 2 points. In the following code adjacent cells are combined.
 718   // If they meet some requirements (~ same curvature, ~ same n) the innermost cell id is added
 719   // to the list of neighbours of the outermost cell. When the cell is added to the neighbours of
 720   // the outermost cell the "level" of the latter is set to the level of the innermost one + 1.
 721   // ( only if $(level of the innermost) + 1 > $(level of the outermost) )
 722   for (int iD = 0; iD < 4; ++iD) {
 723     for (size_t c0 = 0; c0 < fCells[iD].size(); ++c0) {
 724       const int idx = fCells[iD][c0].d1();
 725       for (size_t c1 = tLUT[iD][idx]; idx == fCells[iD + 1][c1].d0(); ++c1) {
 726 #ifdef _TUNING_
 727         fGood = (fLayer[iD].GetClusterSorted(fCells[iD][c0].x())->GetLabel(0) ==
 728                  fLayer[iD + 1].GetClusterSorted(fCells[iD][c0].y())->GetLabel(0) &&
 729                  fLayer[iD + 1].GetClusterSorted(fCells[iD][c0].y())->GetLabel(0) ==
 730                  fLayer[iD + 2].GetClusterSorted(fCells[iD][c0].z())->GetLabel(0) &&
 731                  fLayer[iD + 2].GetClusterSorted(fCells[iD][c0].z())->GetLabel(0) ==
 732                  fLayer[iD + 3].GetClusterSorted(fCells[iD + 1][c1].z())->GetLabel(0) &&
 733                  fLayer[iD].GetClusterSorted(fCells[iD][c0].x())->GetLabel(0) > 0);
 734 #endif
 735         float *n0 = fCells[iD][c0].GetN();
 736         float *n1 = fCells[iD + 1][c1].GetN();
 737         const float dn2 = ((n0[0] - n1[0]) * (n0[0] - n1[0]) + (n0[1] - n1[1]) * (n0[1] - n1[1]) +
 738                            (n0[2] - n1[2]) * (n0[2] - n1[2]));
 739         const float dp = fabs(fCells[iD][c0].GetCurvature() - fCells[iD + 1][c1].GetCurvature());
 740         if (dn2 < fCDN[iD] && dp < fCDP[iD]) {
 741 #ifdef _TUNING_
 742           assert(fCells[iD + 1][c1].Combine(fCells[iD][c0], c0));
 743           if (fGood) {
 744             fGoodCombChi2[iD]->Fill(dp);
 745             fGoodCombN[iD]->Fill(dn2);
 746             cellsCombiningSuccesses++;
 747           } else {
 748             fFakeCombChi2[iD]->Fill(dp);
 749             fFakeCombN[iD]->Fill(dn2);
 750             cellsWrongCombinations++;
 751           }
 752 #else
 753           fCells[iD + 1][c1].Combine(fCells[iD][c0], c0);
 754 #endif
 755         }
 756       }
 757     }
 758   }
 759 #ifdef _TUNING_
 760   Info("MakeCells","Good doublets: %d",numberOfGoodDoublets);
 761   Info("MakeCells","Number of doublets: %d",totalNumberOfDoublets);
 762   Info("MakeCells","Good cells: %d",numberOfGoodCells);
 763   Info("MakeCells","Number of cells: %d",totalNumberOfCells);
 764   Info("MakeCells","Cells combining successes: %d",cellsCombiningSuccesses);
 765   Info("MakeCells","Cells wrong combinations: %d",cellsWrongCombinations);
 766 #endif
 767 }
 768
 769 //__________________________________________________________________________________________________
 770 Int_t AliITSUCATracker::PropagateBack(AliESDEvent * event)
 771 {
 772
 773   Int_t n=event->GetNumberOfTracks();
 774   Int_t ntrk=0;
 775   Int_t ngood=0;
 776   for (Int_t i=0; i<n; i++) {
 777     AliESDtrack *esdTrack=event->GetTrack(i);
 778
 779     if ((esdTrack->GetStatus()&AliESDtrack::kITSin)==0) continue;
 780     if (esdTrack->IsOn(AliESDtrack::kTPCin)) continue; //skip a TPC+ITS track
 781
 782     AliITSUTrackCooked track(*esdTrack);
 783     AliITSUTrackCooked temp(*esdTrack);
 784
 785     temp.ResetCovariance(10.);
 786     temp.ResetClusters();
 787
 788     if (RefitAt(40., &temp, &track)) {
 789
 790       CookLabel(&temp, 0.); //For comparison only
 791       Int_t label = temp.GetLabel();
 792       if (label > 0) ngood++;
 793
 794       esdTrack->UpdateTrackParams(&temp,AliESDtrack::kITSout);
 795       ntrk++;
 796     }
 797   }
 798
 799   Info("PropagateBack","Back propagated tracks: %d",ntrk);
 800   if (ntrk)
 801     Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk);
 802
 803   if (!fSAonly) return AliITSUTrackerGlo::PropagateBack(event);
 804   return 0;
 805 }
 806
 807 //__________________________________________________________________________________________________
 808 Bool_t AliITSUCATracker::RefitAt(Double_t xx, AliITSUTrackCooked *t, const AliITSUTrackCooked *c)
 809 {
 810   // This function refits the track "t" at the position "x" using
 811   // the clusters from "c"
 812
 813   const int nLayers = 7;
 814   Int_t index[nLayers];
 815   Int_t k;
 816   for (k = 0; k < nLayers; k++) index[k] = -1;
 817   Int_t nc = c->GetNumberOfClusters();
 818   for (k = 0; k < nc; k++) {
 819     Int_t idx = c->GetClusterIndex(k), nl = (idx&0xf0000000)>>28;
 820     index[nl] = idx;
 821   }
 822
 823   Int_t from, to, step;
 824   if (xx > t->GetX()) {
 825     from = 0;
 826     to = nLayers;
 827     step = +1;
 828   } else {
 829     from = nLayers - 1;
 830     to = -1;
 831     step = -1;
 832   }
 833
 834   for (Int_t i = from; i != to; i += step) {
 835     Int_t idx = index[i];
 836     if (idx >= 0) {
 837       const AliCluster *cl = GetCluster(idx);
 838       Float_t xr,ar;
 839       cl->GetXAlphaRefPlane(xr, ar);
 840       if (!t->Propagate(Double_t(ar), Double_t(xr), GetBz())) {
 841         return kFALSE;
 842       }
 843       Double_t chi2 = t->GetPredictedChi2(cl);
 844       //      if (chi2 < 100)
 845       t->Update(cl, chi2, idx);
 846     } else {
 847       Double_t r = fgkR[i];
 848       Double_t phi,z;
 849       if (!t->GetPhiZat(r,phi,z)) {
 850         return kFALSE;
 851       }
 852       if (!t->Propagate(phi, r, GetBz())) {
 853         return kFALSE;
 854       }
 855     }
 856     Double_t xx0 = (i > 2) ? 0.008 : 0.003;  // Rough layer thickness
 857     Double_t x0  = 9.36; // Radiation length of Si [cm]
 858     Double_t rho = 2.33; // Density of Si [g/cm^3]
 859     Double_t mass = t->GetMass();
 860     t->CorrectForMeanMaterial(xx0, - step * xx0 * x0 * rho, mass, kTRUE);
 861   }
 862
 863   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
 864   return kTRUE;
 865 }
 866
 867 //__________________________________________________________________________________________________
 868 Int_t AliITSUCATracker::RefitInward(AliESDEvent * event)
 869 {
 870   // Some final refit, after the outliers get removed by the smoother ?
 871   // The clusters must be loaded
 872
 873   Int_t n = event->GetNumberOfTracks();
 874   Int_t ntrk = 0;
 875   Int_t ngood = 0;
 876   for (Int_t i = 0; i < n; i++) {
 877     AliESDtrack *esdTrack = event->GetTrack(i);
 878
 879     if ((esdTrack->GetStatus() & AliESDtrack::kITSout) == 0) continue;
 880     if (esdTrack->IsOn(AliESDtrack::kTPCin)) continue; //skip a TPC+ITS track
 881     AliITSUTrackCooked track(*esdTrack);
 882     AliITSUTrackCooked temp(*esdTrack);
 883
 884     temp.ResetCovariance(10.);
 885     temp.ResetClusters();
 886
 887     if (!RefitAt(2.1, &temp, &track)) continue;
 888     //Cross the beam pipe
 889     if (!temp.PropagateTo(1.8, 2.27e-3, 35.28 * 1.848)) continue;
 890
 891     CookLabel(&temp, 0.); //For comparison only
 892     Int_t label = temp.GetLabel();
 893     if (label > 0) ngood++;
 894
 895     esdTrack->UpdateTrackParams(&temp,AliESDtrack::kITSrefit);
 896     ntrk++;
 897   }
 898
 899   Info("RefitInward","Refitted tracks: %d",ntrk);
 900   if (ntrk)
 901     Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood)/ntrk);
 902
 903   if (!fSAonly) return AliITSUTrackerGlo::RefitInward(event);
 904   return 0;
 905 }
 906
 907 //__________________________________________________________________________________________________
 908 void AliITSUCATracker::SetCuts(int it)
 909 {
 910   switch (it) {
 911     case 0:
 912       fCPhi = fPhiCut;
 913       fCDTanL = kDoublTanL;
 914       fCDPhi = kDoublPhi;
 915       fCZ = fZCut;
 916       for (int i = 0; i < 5; ++i) {
 917         fCDCAxy[i] = kmaxDCAxy[i];
 918         fCDCAz[i] = kmaxDCAz[i];
 919       }
 920       for (int i = 0; i < 4; ++i) {
 921         fCDN[i] = kmaxDN[i];
 922         fCDP[i] = kmaxDP[i];
 923       }
 924       for (int i = 0; i < 6; ++i) {
 925         fCDZ[i] = kmaxDZ[i];
 926       }
 927       break;
 928
 929     default:
 930       fCPhi = 3.f * fPhiCut;
 931       fCDTanL = kDoublTanL1;
 932       fCDPhi = kDoublPhi1;
 933       fCZ = fZCut;
 934       for (int i = 0; i < 5; ++i) {
 935         fCDCAxy[i] = kmaxDCAxy1[i];
 936         fCDCAz[i] = kmaxDCAz1[i];
 937       }
 938       for (int i = 0; i < 4; ++i) {
 939         fCDN[i] = kmaxDN1[i];
 940         fCDP[i] = kmaxDP1[i];
 941       }
 942       for (int i = 0; i < 6; ++i) {
 943         fCDZ[i] = kmaxDZ1[i];
 944       }
 945
 946       break;
 947   }
 948 }
 949
 950 //__________________________________________________________________________________________________
 951 void AliITSUCATracker::UnloadClusters()
 952 {
 953   // This function unloads ITSU clusters from the memory
 954   for (int i = 0;i < 7;++i) {
 955     fLayer[i].Clear();
 956     fUsedClusters[i].clear();
 957   }
 958   for (int i = 0; i < 6; ++i) {
 959     fDoublets[i].clear();
 960   }
 961   for (int i = 0; i < 5; ++i) {
 962     fCells[i].clear();
 963   }
 964   for (int i = 0; i < 4; ++i)
 965   {
 966     fCandidates[i]->Clear("C");
 967   }
 968   AliITSUTrackerGlo::UnloadClusters();
 969 }
 970
 971 #ifdef _TUNING_
 972 //__________________________________________________________________________________________________
 973 void AliITSUCATracker::ResetHistos()
 974 {
 975   for (int i = 0; i < 6; ++i) {
 976     fGDZ[i]->Reset();
 977     fGDXY[i]->Reset();
 978     fFDZ[i]->Reset();
 979     fFDXY[i]->Reset();
 980   }
 981   for (int i = 0; i < 5; ++i) {
 982     fGoodCombChi2[i]->Reset();
 983     fFakeCombChi2[i]->Reset();
 984     fGDCAZ[i]->Reset();
 985     fGDCAXY[i]->Reset();
 986     fFDCAZ[i]->Reset();
 987     fFDCAXY[i]->Reset();
 988   }
 989   for (int i = 0; i < 4; ++i) {
 990     fGoodCombN[i]->Reset();
 991     fFakeCombN[i]->Reset();
 992   }
 993   fTan->Reset();
 994   fTanF->Reset();
 995   fPhi->Reset();
 996   fPhiF->Reset();
 997   fNEntries->Reset();
 998 }
 999 #endif