30280940a34ea605491836856eb977a3853ea9f4
[u/mrichter/AliRoot.git] / 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