Integrating the Cooked Matrix tracker into the commom reconstruction framework
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUTrackerCooked.cxx
1 //-------------------------------------------------------------------------
2 //               Implementation of the ITS tracker class
3 //    The pattern recongintion based on the "cooked covariance" approach
4 //-------------------------------------------------------------------------
5
6 #include <TTree.h>
7 #include <TClonesArray.h>
8
9 #include "AliLog.h"
10 #include "AliESDEvent.h"
11 #include "AliITSUClusterPix.h"
12 #include "AliITSUGeomTGeo.h"
13 #include "AliITSUTrackerCooked.h"
14 #include "AliITSUTrackCooked.h" 
15 #include "AliITSUReconstructor.h" 
16
17 ClassImp(AliITSUTrackerCooked)
18
19 //************************************************
20 // Constants hardcoded for the moment:
21 //************************************************
22 // seed "windows" in z and phi: MakeSeeds
23 const Double_t kzWin=0.33;
24 const Double_t kminPt=0.05;
25 // Maximal accepted impact parameters for the seeds 
26 const Double_t kmaxDCAxy=3.;
27 const Double_t kmaxDCAz= 3.;
28 // Layers for the seeding
29 const Int_t kSeedingLayer1=6, kSeedingLayer2=4, kSeedingLayer3=5;
30 // Space point resolution
31 const Double_t kSigma2=0.0005*0.0005;
32 // Max accepted chi2
33 const Double_t kmaxChi2PerCluster=20.;
34 const Double_t kmaxChi2PerTrack=30.;
35 // Tracking "road" from layer to layer
36 const Double_t kRoadY=0.7;
37 const Double_t kRoadZ=0.7;
38 // Minimal number of attached clusters
39 const Int_t kminNumberOfClusters=4;
40
41 //************************************************
42 // TODO:
43 //************************************************
44 // Seeding:
45 // Precalculate cylidnrical (r,phi) for the clusters;
46 // use exact r's for the clusters
47
48
49 AliITSUTrackerCooked::AliITSUlayer
50               AliITSUTrackerCooked::fgLayers[AliITSUTrackerCooked::kNLayers];
51
52 AliITSUTrackerCooked::AliITSUTrackerCooked(AliITSUReconstructor *rec): 
53 AliITSUTrackerGlo(rec),
54 fSeeds(0),
55 fI(kNLayers-1),
56 fBestTrack(0), 
57 fTrackToFollow(0),
58 fSAonly(kTRUE) 
59 {
60   //--------------------------------------------------------------------
61   // This default constructor needs to be provided
62   //--------------------------------------------------------------------
63   const Double_t 
64   klRadius[7]={2.34, 3.15, 3.93, 19.61, 24.55, 34.39, 39.34}; //tdr6
65
66   AliITSUGeomTGeo *gm  = new AliITSUGeomTGeo(kTRUE,kTRUE);
67   AliITSUClusterPix::SetGeom(gm);
68
69   for (Int_t i=0; i<kNLayers; i++) fgLayers[i].SetR(klRadius[i]);
70
71   // Some default primary vertex
72   Double_t xyz[]={0.,0.,0.};
73   Double_t ers[]={2.,2.,2.};
74
75   SetVertex(xyz,ers);
76
77 }
78
79 void AliITSUTrackerCooked::ResetTrackToFollow(const AliITSUTrackCooked &t) {
80   //--------------------------------------------------------------------
81   // Prepare to follow a new track seed
82   //--------------------------------------------------------------------
83      delete fTrackToFollow;
84      fTrackToFollow = new AliITSUTrackCooked(t);
85 }
86   
87 void AliITSUTrackerCooked::ResetBestTrack() {
88   //--------------------------------------------------------------------
89   // Replace the best track branch
90   //--------------------------------------------------------------------
91      delete fBestTrack;
92      fBestTrack = new AliITSUTrackCooked(*fTrackToFollow);
93 }
94   
95 AliITSUTrackerCooked::~AliITSUTrackerCooked() 
96 {
97   //--------------------------------------------------------------------
98   // Virtual destructor
99   //--------------------------------------------------------------------
100
101   if (fSeeds) fSeeds->Delete(); delete fSeeds; 
102   delete fBestTrack;
103   delete fTrackToFollow;
104
105 }
106
107 static Double_t 
108 f1(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
109 {
110     //-----------------------------------------------------------------
111     // Initial approximation of the track curvature
112     //-----------------------------------------------------------------
113     Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
114     Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
115                     (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
116     Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
117                     (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
118     
119     Double_t xr=TMath::Abs(d/(d*x1-a)), yr=TMath::Abs(d/(d*y1-b));
120     
121     Double_t crv=xr*yr/sqrt(xr*xr+yr*yr);
122     if (d>0) crv=-crv;
123
124     return crv;
125 }
126
127 static Double_t 
128 f2(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
129 {
130     //-----------------------------------------------------------------
131     // Initial approximation of the x-coordinate of the center of curvature 
132     //-----------------------------------------------------------------
133
134   Double_t k1=(y2-y1)/(x2-x1), k2=(y3-y2)/(x3-x2);
135   Double_t x0=0.5*(k1*k2*(y1-y3) + k2*(x1+x2) - k1*(x2+x3))/(k2-k1);
136
137   return x0;
138 }
139
140 static Double_t 
141 f3(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t z1, Double_t z2)
142 {
143     //-----------------------------------------------------------------
144     // Initial approximation of the tangent of the track dip angle
145     //-----------------------------------------------------------------
146     return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
147 }
148
149 Bool_t AliITSUTrackerCooked::
150 AddCookedSeed(const Float_t r1[3], Int_t l1, Int_t i1, 
151               const Float_t r2[3], Int_t l2, Int_t i2,
152               const AliCluster *c3,Int_t l3, Int_t i3) 
153 {
154     //--------------------------------------------------------------------
155     // This is the main cooking function.
156     // Creates seed parameters out of provided clusters.
157     //--------------------------------------------------------------------
158     Float_t x,a;
159     if (!c3->GetXAlphaRefPlane(x,a)) return kFALSE;
160
161     Double_t ca=TMath::Cos(a), sa=TMath::Sin(a);
162     Double_t x1 = r1[0]*ca + r1[1]*sa,
163              y1 =-r1[0]*sa + r1[1]*ca, z1 = r1[2];
164     Double_t x2 = r2[0]*ca + r2[1]*sa,
165              y2 =-r2[0]*sa + r2[1]*ca, z2 = r2[2];
166     Double_t x3 = x,  y3 = c3->GetY(), z3 = c3->GetZ();
167
168     Double_t par[5];
169     par[0]=y3;
170     par[1]=z3;
171     Double_t crv=f1(x1, y1, x2, y2, x3, y3); //curvature
172     Double_t x0 =f2(x1, y1, x2, y2, x3, y3); //x-coordinate of the center
173     Double_t tgl12=f3(x1, y1, x2, y2, z1, z2);
174     Double_t tgl23=f3(x2, y2, x3, y3, z2, z3);
175
176     Double_t sf=crv*(x-x0);
177     if (TMath::Abs(sf) >= kAlmost1) return kFALSE;
178     par[2]=sf;
179
180     par[3]=0.5*(tgl12 + tgl23);
181     Double_t bz=GetBz();
182     par[4]=(TMath::Abs(bz) < kAlmost0Field) ? kAlmost0 : crv/(bz*kB2C);
183
184     Double_t cov[15];
185     /*
186     for (Int_t i=0; i<15; i++) cov[i]=0.;
187     cov[0] =kSigma2*10;
188     cov[2] =kSigma2*10;
189     cov[5] =0.007*0.007*10;   //FIXME all these lines
190     cov[9] =0.007*0.007*10;
191     cov[14]=0.1*0.1*10;
192     */
193     const Double_t dlt=0.0005;
194     Double_t
195     fy=1./(fgLayers[kSeedingLayer3].GetR() - fgLayers[kSeedingLayer2].GetR());
196     Double_t tz=fy;
197     Double_t cy=(f1(x1, y1, x2, y2+dlt, x3, y3) - crv)/dlt/bz/kB2C;
198     cy*=20; //FIXME: MS contribution to the cov[14]
199     Double_t s2=kSigma2;
200
201     cov[0]=s2;
202     cov[1]=0.;     cov[2]=s2;
203     cov[3]=s2*fy;  cov[4]=0.;    cov[5]=s2*fy*fy;
204     cov[6]=0.;     cov[7]=s2*tz; cov[8]=0.;        cov[9]=s2*tz*tz;
205     cov[10]=s2*cy; cov[11]=0.;   cov[12]=s2*fy*cy; cov[13]=0.; cov[14]=s2*cy*cy;
206
207     AliITSUTrackCooked *seed=new AliITSUTrackCooked();
208     seed->Set(Double_t(x), Double_t(a), par, cov);
209
210     Float_t dz[2]; 
211     seed->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
212     if (TMath::Abs(dz[0]) > kmaxDCAxy) {delete seed; return kFALSE;} 
213     if (TMath::Abs(dz[1]) > kmaxDCAz ) {delete seed; return kFALSE;} 
214
215     Double_t xx0 = 0.008; // Rough layer thickness
216     Double_t radl= 9.36;  // Radiation length of Si [cm]
217     Double_t rho = 2.33;  // Density of Si [g/cm^3] 
218     Double_t mass= 0.139;// Pion
219     if (!seed->CorrectForMeanMaterial(xx0, xx0*radl*rho, mass, kTRUE)) {
220        delete seed; return kFALSE;
221     }
222
223     seed->SetClusterIndex(l1,i1);
224     seed->SetClusterIndex(l2,i2);
225     seed->SetClusterIndex(l3,i3);
226
227     fSeeds->AddLast(seed);
228
229     return kTRUE;
230 }
231
232 Int_t AliITSUTrackerCooked::MakeSeeds() {
233   //--------------------------------------------------------------------
234   // This is the main pattern recongition function.
235   // Creates seeds out of two clusters and another point.
236   //--------------------------------------------------------------------
237    if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
238    fSeeds=new TObjArray(77777);
239
240    const Double_t zv=GetZ();
241
242    AliITSUlayer &layer1=fgLayers[kSeedingLayer1];
243    AliITSUlayer &layer2=fgLayers[kSeedingLayer2];
244    AliITSUlayer &layer3=fgLayers[kSeedingLayer3];
245    Double_t r1=layer1.GetR();
246    Double_t r2=layer2.GetR();
247    Double_t r3=layer3.GetR();
248
249    const Double_t maxC  = TMath::Abs(GetBz()*kB2C/kminPt);
250    const Double_t kpWin = TMath::ASin(0.5*maxC*r1) - TMath::ASin(0.5*maxC*r2);
251
252    Int_t nClusters1=layer1.GetNumberOfClusters();
253    Int_t nClusters2=layer2.GetNumberOfClusters();
254    Int_t nClusters3=layer3.GetNumberOfClusters();
255    for (Int_t n1=0; n1<nClusters1; n1++) {
256      AliCluster *c1=layer1.GetCluster(n1);
257      //
258      //Int_t lab=c1->GetLabel(0);
259      //
260      Double_t z1=c1->GetZ();
261      Float_t xyz1[3]; c1->GetGlobalXYZ(xyz1);
262      Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
263      Double_t zr2=zv + r2/r1*(z1-zv);
264      Int_t start2=layer2.FindClusterIndex(zr2-kzWin);
265      for (Int_t n2=start2; n2<nClusters2; n2++) {
266          AliCluster *c2=layer2.GetCluster(n2);
267          //
268          //if (c2->GetLabel(0)!=lab) continue;
269          //
270          Double_t z2=c2->GetZ();
271          if (z2 > (zr2+kzWin)) break;  //check in Z
272
273          Float_t xyz2[3]; c2->GetGlobalXYZ(xyz2);
274          Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
275          if (TMath::Abs(phi2-phi1) > kpWin) continue;  //check in Phi
276
277          Double_t zr3=z1 + (r3-r1)/(r2-r1)*(z2-z1);
278          Double_t crv=f1(xyz1[0], xyz1[1], xyz2[0], xyz2[1], GetX(), GetY());
279          Double_t phir3 = phi1 + 0.5*crv*(r3 - r1); 
280
281          Int_t start3=layer3.FindClusterIndex(zr3-kzWin/2);
282          for (Int_t n3=start3; n3<nClusters3; n3++) {
283              AliCluster *c3=layer3.GetCluster(n3);
284              //
285              //if (c3->GetLabel(0)!=lab) continue;
286              //
287              Double_t z3=c3->GetZ();
288              if (z3 > (zr3+kzWin/2)) break;  //check in Z
289
290              Float_t xyz3[3]; c3->GetGlobalXYZ(xyz3);
291              Double_t phi3=TMath::ATan2(xyz3[1],xyz3[0]);
292              if (TMath::Abs(phir3-phi3) > kpWin/100) continue;  //check in Phi
293
294              AliITSUClusterPix cc(*((AliITSUClusterPix*)c2));
295              cc.GoToFrameTrk();
296              AddCookedSeed(xyz1, kSeedingLayer1, n1,
297                            xyz3, kSeedingLayer3, n3, 
298                            &cc,  kSeedingLayer2, n2);
299
300          }
301      }
302    }
303
304    for (Int_t n1=0; n1<nClusters1; n1++) {
305      AliCluster *c1=layer1.GetCluster(n1);
306      ((AliITSUClusterPix*)c1)->GoToFrameTrk();
307    }
308    for (Int_t n2=0; n2<nClusters2; n2++) {
309      AliCluster *c2=layer2.GetCluster(n2);
310      ((AliITSUClusterPix*)c2)->GoToFrameTrk();
311    }
312    for (Int_t n3=0; n3<nClusters3; n3++) {
313      AliCluster *c3=layer3.GetCluster(n3);
314      ((AliITSUClusterPix*)c3)->GoToFrameTrk();
315    }
316
317    fSeeds->Sort();
318    return fSeeds->GetEntriesFast();
319 }
320
321 Int_t AliITSUTrackerCooked::Clusters2Tracks(AliESDEvent *event) {
322   //--------------------------------------------------------------------
323   // This is the main tracking function
324   // The clusters must already be loaded
325   //--------------------------------------------------------------------
326
327   if (!fSAonly) AliITSUTrackerGlo::Clusters2Tracks(event);
328
329   Int_t nSeeds=MakeSeeds();
330
331   // Possibly, icrement the seeds with additional clusters (Kalman)
332
333   // Possibly, (re)fit the found tracks 
334
335   Int_t ngood=0;
336   for (Int_t s=0; s<nSeeds; s++) {
337       const AliITSUTrackCooked *track=(AliITSUTrackCooked*)fSeeds->At(s);
338       ResetTrackToFollow(*track);
339       ResetBestTrack();
340       fI=kSeedingLayer2;
341       fgLayers[fI].ResetTrack(*track);
342
343       for (FollowProlongation(); fI<kSeedingLayer2; fI++) {
344           while (TakeNextProlongation()) FollowProlongation();
345       }
346
347       if (fBestTrack->GetNumberOfClusters() < kminNumberOfClusters) continue;
348
349       CookLabel(fBestTrack,0.); //For comparison only
350       Int_t label=fBestTrack->GetLabel();
351       if (label>0) ngood++;
352
353       AliESDtrack iotrack;
354       iotrack.UpdateTrackParams(fBestTrack,AliESDtrack::kITSin);
355       iotrack.SetLabel(label);
356       event->AddTrack(&iotrack);
357       UseClusters(fBestTrack);
358   }
359
360   Info("Clusters2Tracks","Seeds: %d",nSeeds);
361   if (nSeeds)
362   Info("Clusters2Tracks","Good tracks/seeds: %f",Float_t(ngood)/nSeeds);
363
364   if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
365   fSeeds=0;
366     
367   return 0;
368 }
369
370 void AliITSUTrackerCooked::FollowProlongation() {
371   //--------------------------------------------------------------------
372   // Push this track tree branch towards the primary vertex
373   //--------------------------------------------------------------------
374   while (fI) {
375     fI--;
376     AliITSUlayer &layer = fgLayers[fI]; //fI is the number of the next layer
377     Double_t r=layer.GetR();
378
379     //Find intersection (fTrackToFollow is still at the previous layer)
380     Double_t phi,z;  
381     if (!fTrackToFollow->GetPhiZat(r,phi,z)) {
382       //Warning("FollowProlongation","failed to estimate track !\n");
383       return;
384     }
385
386     //if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) return;
387     Double_t zMin = z - kRoadZ; 
388     Double_t zMax = z + kRoadZ;
389     Double_t phiMin = phi - kRoadY/r;
390     Double_t phiMax = phi + kRoadY/r;
391     if (layer.SelectClusters(zMin, zMax, phiMin, phiMax)==0) return;  
392
393     if (!TakeNextProlongation()) return;
394
395   } 
396
397   //deal with the best track
398   Int_t ncl=fTrackToFollow->GetNumberOfClusters();
399   Int_t nclb=fBestTrack->GetNumberOfClusters();
400   if (ncl >= nclb) {
401      Double_t chi2=fTrackToFollow->GetChi2();
402      if (chi2 < kmaxChi2PerTrack) {        
403         if (ncl > nclb || chi2 < fBestTrack->GetChi2()) {
404            ResetBestTrack();
405         }
406      }
407   }
408
409 }
410
411 Int_t AliITSUTrackerCooked::TakeNextProlongation() {
412   //--------------------------------------------------------------------
413   // Switch to the next track tree branch
414   //--------------------------------------------------------------------
415   AliITSUlayer &layer=fgLayers[fI];
416
417   const AliCluster *c=0; Int_t ci=-1;
418   const AliCluster *cc=0; Int_t cci=-1;
419   UShort_t volId=-1;
420   Double_t x=0., alpha=0.;
421   Double_t z=0., dz=0., y=0., dy=0., chi2=0.; 
422   while ((c=layer.GetNextCluster(ci))!=0) {
423     if (c->IsClusterUsed()) continue;
424     Int_t id=c->GetVolumeId();
425     if (id != volId) {
426        volId=id;
427        Float_t xr,ar; c->GetXAlphaRefPlane(xr, ar);
428        x=xr; alpha=ar;
429        const AliITSUTrackCooked *t = fgLayers[fI+1].GetTrack();
430        ResetTrackToFollow(*t);
431        if (!fTrackToFollow->Propagate(alpha, x, GetBz())) {
432          //Warning("TakeNextProlongation","propagation failed !\n");
433           continue;
434        }
435        dz=7*TMath::Sqrt(fTrackToFollow->GetSigmaZ2() + kSigma2);
436        dy=7*TMath::Sqrt(fTrackToFollow->GetSigmaY2() + kSigma2);
437        z=fTrackToFollow->GetZ();
438        y=fTrackToFollow->GetY();
439     }
440
441     //if (TMath::Abs(fTrackToFollow.GetZ()-GetZ())>layer.GetR()+dz) continue;
442
443     if (TMath::Abs(z - c->GetZ()) > dz) continue;
444     if (TMath::Abs(y - c->GetY()) > dy) continue;
445
446     Double_t ch2=fTrackToFollow->GetPredictedChi2(c); 
447     if (ch2 > kmaxChi2PerCluster) continue;
448     chi2=ch2;
449     cc=c; cci=ci;
450     break;
451   }
452
453   if (!cc) return 0;
454
455   if (!fTrackToFollow->Update(cc,chi2,(fI<<28)+cci)) {
456      //Warning("TakeNextProlongation","filtering failed !\n");
457      return 0;
458   }
459   Double_t xx0 = (fI > 2) ? 0.008 : 0.003;  // Rough layer thickness
460   Double_t x0  = 9.36; // Radiation length of Si [cm]
461   Double_t rho = 2.33; // Density of Si [g/cm^3] 
462   Double_t mass = fTrackToFollow->GetMass();
463   fTrackToFollow->CorrectForMeanMaterial(xx0, xx0*x0*rho, mass, kTRUE);
464   layer.ResetTrack(*fTrackToFollow); 
465
466   return 1;
467 }
468
469 Int_t AliITSUTrackerCooked::PropagateBack(AliESDEvent *event) {
470   //--------------------------------------------------------------------
471   // Here, we implement the Kalman smoother ?
472   // The clusters must already be loaded
473   //--------------------------------------------------------------------
474   Int_t n=event->GetNumberOfTracks();
475   Int_t ntrk=0;
476   Int_t ngood=0;
477   for (Int_t i=0; i<n; i++) {
478       AliESDtrack *esdTrack=event->GetTrack(i);
479
480       if (!esdTrack->IsOn(AliESDtrack::kITSin)) continue;
481       if ( esdTrack->IsOn(AliESDtrack::kTPCin)) continue;//skip a TPC+ITS track
482
483       AliITSUTrackCooked track(*esdTrack);
484
485       ResetTrackToFollow(track);
486
487       fTrackToFollow->ResetCovariance(10.); fTrackToFollow->ResetClusters();
488       if (RefitAt(40., fTrackToFollow, &track)) {
489
490          CookLabel(fTrackToFollow, 0.); //For comparison only
491          Int_t label=fTrackToFollow->GetLabel();
492          if (label>0) ngood++;
493
494          esdTrack->UpdateTrackParams(fTrackToFollow,AliESDtrack::kITSout);
495          //UseClusters(fTrackToFollow);
496          ntrk++;
497       }
498   }
499
500   Info("PropagateBack","Back propagated tracks: %d",ntrk);
501   if (ntrk)
502   Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk);
503   
504   if (!fSAonly) AliITSUTrackerGlo::PropagateBack(event);
505   
506   return 0;
507 }
508
509 Bool_t AliITSUTrackerCooked::
510 RefitAt(Double_t xx, AliITSUTrackCooked *t, const AliITSUTrackCooked *c) {
511   //--------------------------------------------------------------------
512   // This function refits the track "t" at the position "x" using
513   // the clusters from "c"
514   //--------------------------------------------------------------------
515   Int_t index[kNLayers];
516   Int_t k;
517   for (k=0; k<kNLayers; k++) index[k]=-1;
518   Int_t nc=c->GetNumberOfClusters();
519   for (k=0; k<nc; k++) {
520     Int_t idx=c->GetClusterIndex(k), nl=(idx&0xf0000000)>>28;
521     index[nl]=idx;
522   }
523
524   Int_t from, to, step;
525   if (xx > t->GetX()) {
526       from=0; to=kNLayers;
527       step=+1;
528   } else {
529       from=kNLayers-1; to=-1;
530       step=-1;
531   }
532
533   for (Int_t i=from; i != to; i += step) {
534      Int_t idx=index[i];
535      if (idx>=0) {
536         const AliCluster *cl=GetCluster(idx);
537         Float_t xr,ar; cl->GetXAlphaRefPlane(xr, ar);
538         if (!t->Propagate(Double_t(ar), Double_t(xr), GetBz())) {
539            //Warning("RefitAt","propagation failed !\n");
540            return kFALSE;
541         }
542         Double_t chi2=t->GetPredictedChi2(cl);
543         if (chi2 < kmaxChi2PerCluster) t->Update(cl, chi2, idx);
544      } else {
545         Double_t r=fgLayers[i].GetR();
546         Double_t phi,z;
547         if (!t->GetPhiZat(r,phi,z)) {
548            //Warning("RefitAt","failed to estimate track !\n");
549            return kFALSE;
550         }
551         if (!t->Propagate(phi, r, GetBz())) {
552            //Warning("RefitAt","propagation failed !\n");
553            return kFALSE;
554         }
555      }
556      Double_t xx0 = (i > 2) ? 0.008 : 0.003;  // Rough layer thickness
557      Double_t x0  = 9.36; // Radiation length of Si [cm]
558      Double_t rho = 2.33; // Density of Si [g/cm^3]
559      Double_t mass = t->GetMass();
560      t->CorrectForMeanMaterial(xx0, -step*xx0*x0*rho, mass, kTRUE);
561   }
562
563   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
564   return kTRUE;
565 }
566
567 Int_t AliITSUTrackerCooked::RefitInward(AliESDEvent *event) {
568   //--------------------------------------------------------------------
569   // Some final refit, after the outliers get removed by the smoother ?  
570   // The clusters must be loaded
571   //--------------------------------------------------------------------
572   Int_t n=event->GetNumberOfTracks();
573   Int_t ntrk=0;
574   Int_t ngood=0;
575   for (Int_t i=0; i<n; i++) {
576       AliESDtrack *esdTrack=event->GetTrack(i);
577
578       if (!esdTrack->IsOn(AliESDtrack::kITSout)) continue;
579       if ( esdTrack->IsOn(AliESDtrack::kTPCin)) continue;//skip a TPC+ITS track
580
581       AliITSUTrackCooked track(*esdTrack);
582       ResetTrackToFollow(track);
583
584       fTrackToFollow->ResetCovariance(10.); fTrackToFollow->ResetClusters();
585       if (!RefitAt(2.1, fTrackToFollow, &track)) continue;
586       //Cross the beam pipe
587       if (!fTrackToFollow->PropagateTo(1.8, 2.27e-3, 35.28*1.848)) continue;
588
589       CookLabel(fTrackToFollow, 0.); //For comparison only
590       Int_t label=fTrackToFollow->GetLabel();
591       if (label>0) ngood++;
592
593       esdTrack->UpdateTrackParams(fTrackToFollow,AliESDtrack::kITSrefit);
594       //esdTrack->RelateToVertex(event->GetVertex(),GetBz(),33.);
595       //UseClusters(fTrackToFollow);
596       ntrk++;
597   }
598
599   Info("RefitInward","Refitted tracks: %d",ntrk);
600   if (ntrk)
601   Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood)/ntrk);
602     
603   if (!fSAonly) AliITSUTrackerGlo::RefitInward(event);
604
605   return 0;
606 }
607
608 Int_t AliITSUTrackerCooked::LoadClusters(TTree *cTree) {
609   //--------------------------------------------------------------------
610   // This function reads the ITSU clusters from the tree,
611   // sort them, distribute over the internal tracker arrays, etc
612   //--------------------------------------------------------------------
613   if (!cTree) {
614      AliFatal("No cluster tree !");
615      return 1;
616   }
617
618   AliITSUTrackerGlo::LoadClusters(cTree);
619
620   for (Int_t i=0; i<kNLayers; i++) {
621       TClonesArray *clusters=fReconstructor->GetClusters(i);
622       switch (i) {
623       case kSeedingLayer1: 
624       case kSeedingLayer2: 
625       case kSeedingLayer3: 
626          fgLayers[i].InsertClusters(clusters,kTRUE,fSAonly);
627          break;
628       default:
629          fgLayers[i].InsertClusters(clusters,kFALSE,fSAonly);
630          break;
631       }
632   }
633
634   return 0;
635 }
636
637 void AliITSUTrackerCooked::UnloadClusters() {
638   //--------------------------------------------------------------------
639   // This function unloads ITSU clusters from the RAM
640   //--------------------------------------------------------------------
641   AliITSUTrackerGlo::UnloadClusters();
642   for (Int_t i=0; i<kNLayers; i++) fgLayers[i].DeleteClusters();
643 }
644
645 AliCluster *AliITSUTrackerCooked::GetCluster(Int_t index) const {
646   //--------------------------------------------------------------------
647   //       Return pointer to a given cluster
648   //--------------------------------------------------------------------
649     Int_t l=(index & 0xf0000000) >> 28;
650     Int_t c=(index & 0x0fffffff) >> 00;
651     return fgLayers[l].GetCluster(c);
652 }
653
654 AliITSUTrackerCooked::AliITSUlayer::AliITSUlayer():
655   fR(0),
656   fN(0),
657   fNsel(0),
658   fTrack(0) 
659 {
660   //--------------------------------------------------------------------
661   // This default constructor needs to be provided
662   //--------------------------------------------------------------------
663   for (Int_t i=0; i<kMaxClusterPerLayer; i++) fClusters[i]=0;
664   for (Int_t i=0; i<kMaxSelected; i++) fIndex[i]=-1;
665 }
666
667 AliITSUTrackerCooked::AliITSUlayer::~AliITSUlayer()
668 {
669   //--------------------------------------------------------------------
670   // Simple destructor
671   //--------------------------------------------------------------------
672   delete fTrack;
673 }
674
675 void 
676 AliITSUTrackerCooked::AliITSUlayer::ResetTrack(const AliITSUTrackCooked &t) {
677   //--------------------------------------------------------------------
678   // Replace the track estimate at this layer
679   //--------------------------------------------------------------------
680    delete fTrack;
681    fTrack=new AliITSUTrackCooked(t);
682 }
683
684 void AliITSUTrackerCooked::AliITSUlayer::
685 InsertClusters(TClonesArray *clusters, Bool_t seedingLayer, Bool_t saOnly)
686 {
687   //--------------------------------------------------------------------
688   // Load clusters to this layer
689   //--------------------------------------------------------------------
690   Int_t ncl=clusters->GetEntriesFast();
691   Double_t r=0.;
692   for (Int_t i=0; i<ncl; i++) {
693      AliITSUClusterPix *c=(AliITSUClusterPix*)clusters->UncheckedAt(i);
694      if (!saOnly) if (c->IsClusterUsed()) continue;
695      c->GoToFrameGlo();
696      Double_t x=c->GetX(), y=c->GetY();
697      r += TMath::Sqrt(x*x + y*y);
698      if (!seedingLayer) c->GoToFrameTrk();
699      //if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
700      InsertCluster(c);
701   }
702   if (ncl) fR = r/ncl;
703 }
704
705 void AliITSUTrackerCooked::AliITSUlayer::DeleteClusters()
706 {
707   //--------------------------------------------------------------------
708   // Load clusters to this layer
709   //--------------------------------------------------------------------
710   //for (Int_t i=0; i<fN; i++) {delete fClusters[i]; fClusters[i]=0;}
711   fN=0;
712 }
713
714 Int_t 
715 AliITSUTrackerCooked::AliITSUlayer::InsertCluster(AliCluster *c) {
716   //--------------------------------------------------------------------
717   // This function inserts a cluster to this layer in increasing
718   // order of the cluster's fZ
719   //--------------------------------------------------------------------
720   if (fN>=kMaxClusterPerLayer) {
721      ::Error("InsertCluster","Too many clusters !\n");
722      return 1;
723   }
724   if (fN==0) fClusters[0]=c;
725   else {
726      Int_t i=FindClusterIndex(c->GetZ());
727      Int_t k=fN-i;
728      memmove(fClusters+i+1 ,fClusters+i,k*sizeof(AliCluster*));
729      fClusters[i]=c;
730   }
731   fN++;
732   return 0;
733 }
734
735 Int_t 
736 AliITSUTrackerCooked::AliITSUlayer::FindClusterIndex(Double_t z) const {
737   //--------------------------------------------------------------------
738   // This function returns the index of the first 
739   // with its fZ >= "z". 
740   //--------------------------------------------------------------------
741   if (fN==0) return 0;
742
743   Int_t b=0;
744   if (z <= fClusters[b]->GetZ()) return b;
745
746   Int_t e=b+fN-1;
747   if (z > fClusters[e]->GetZ()) return e+1;
748
749   Int_t m=(b+e)/2;
750   for (; b<e; m=(b+e)/2) {
751     if (z > fClusters[m]->GetZ()) b=m+1;
752     else e=m; 
753   }
754   return m;
755 }
756
757 Int_t AliITSUTrackerCooked::AliITSUlayer::
758 SelectClusters(Float_t zMin,Float_t zMax,Float_t phiMin, Float_t phiMax) {
759   //--------------------------------------------------------------------
760   // This function selects clusters within the "road"
761   //--------------------------------------------------------------------
762   UShort_t volId=-1;
763   Float_t x=0., alpha=0.;
764   for (Int_t i=FindClusterIndex(zMin); i<fN; i++) {
765       AliCluster *c=fClusters[i];
766       if (c->GetZ() > zMax) break;
767       if (c->IsClusterUsed()) continue;
768       UShort_t id=c->GetVolumeId();
769       if (id != volId) {
770         volId=id;
771         c->GetXAlphaRefPlane(x,alpha); //FIXME
772       }
773       Double_t cPhi=alpha + c->GetY()/fR;
774       if (cPhi<0.) cPhi+=2.*TMath::Pi();
775       else if (cPhi >= 2.*TMath::Pi()) cPhi-=2.*TMath::Pi();
776       if (cPhi <= phiMin) continue;
777       if (cPhi >  phiMax) continue;
778       fIndex[fNsel++]=i;
779       if (fNsel==kMaxSelected) break;
780   }
781   return fNsel;
782 }
783
784 const AliCluster *AliITSUTrackerCooked::AliITSUlayer::GetNextCluster(Int_t &ci){
785   //--------------------------------------------------------------------
786   // This function returns clusters within the "road" 
787   //--------------------------------------------------------------------
788   AliCluster *c=0;
789   ci=-1;
790   if (fNsel) {
791      fNsel--;
792      ci=fIndex[fNsel]; 
793      c=fClusters[ci];
794   }
795   return c; 
796 }
797