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