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