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=layer1.GetClusterPhi(n1);
 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=layer2.GetClusterPhi(n2);
 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=layer3.GetClusterPhi(n3);
 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
 339       Double_t x=track->GetX();
 340       Double_t y=track->GetY();
 341       Double_t phi=track->GetAlpha() + TMath::ATan2(y,x);
 342       const Float_t pi2 = 2.*TMath::Pi();
 343       if (phi<0.) phi+=pi2;
 344       else if (phi >= pi2) phi-=pi2;
 345       for (Int_t n=0; n<kNLayers-3; n++) {
 346         Double_t z;
 347         track->GetZAt(fgLayers[n].GetR(),GetBz(),z);
 348         fgLayers[n].SelectClusters(phi,kRoadY,z,kRoadZ);
 349       }
 350
 351       ResetTrackToFollow(*track);
 352       ResetBestTrack();
 353       fI=kSeedingLayer2;
 354       fgLayers[fI].ResetTrack(*track);
 355
 356       for (FollowProlongation(); fI<kSeedingLayer2; fI++) {
 357           while (TakeNextProlongation()) FollowProlongation();
 358       }
 359
 360       if (fBestTrack->GetNumberOfClusters() < kminNumberOfClusters) continue;
 361
 362       CookLabel(fBestTrack,0.); //For comparison only
 363       Int_t label=fBestTrack->GetLabel();
 364       if (label>0) ngood++;
 365
 366       AliESDtrack iotrack;
 367       iotrack.UpdateTrackParams(fBestTrack,AliESDtrack::kITSin);
 368       iotrack.SetLabel(label);
 369       event->AddTrack(&iotrack);
 370       UseClusters(fBestTrack);
 371   }
 372
 373   Info("Clusters2Tracks","Seeds: %d",nSeeds);
 374   if (nSeeds)
 375   Info("Clusters2Tracks","Good tracks/seeds: %f",Float_t(ngood)/nSeeds);
 376
 377   if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
 378   fSeeds=0;
 379
 380   return 0;
 381 }
 382
 383 void AliITSUTrackerCooked::FollowProlongation() {
 384   //--------------------------------------------------------------------
 385   // Push this track tree branch towards the primary vertex
 386   //--------------------------------------------------------------------
 387   while (fI) {
 388     fI--;
 389     fgLayers[fI].ResetSelectedClusters();
 390     if (!TakeNextProlongation()) return;
 391   }
 392
 393   //deal with the best track
 394   Int_t ncl=fTrackToFollow->GetNumberOfClusters();
 395   Int_t nclb=fBestTrack->GetNumberOfClusters();
 396   if (ncl >= nclb) {
 397      Double_t chi2=fTrackToFollow->GetChi2();
 398      if (chi2 < kmaxChi2PerTrack) {
 399         if (ncl > nclb || chi2 < fBestTrack->GetChi2()) {
 400            ResetBestTrack();
 401         }
 402      }
 403   }
 404
 405 }
 406
 407 Int_t AliITSUTrackerCooked::TakeNextProlongation() {
 408   //--------------------------------------------------------------------
 409   // Switch to the next track tree branch
 410   //--------------------------------------------------------------------
 411   AliITSUlayer &layer=fgLayers[fI];
 412
 413   const AliCluster *c=0; Int_t ci=-1;
 414   const AliCluster *cc=0; Int_t cci=-1;
 415   UShort_t volId=-1;
 416   Double_t z=0., dz=0., y=0., dy=0., chi2=0.;
 417   while ((c=layer.GetNextCluster(ci))!=0) {
 418     Int_t id=c->GetVolumeId();
 419     if (id != volId) {
 420        volId=id;
 421        const AliITSUTrackCooked *t = fgLayers[fI+1].GetTrack();
 422        ResetTrackToFollow(*t);
 423        Double_t x=layer.GetXRef(ci);
 424        Double_t alpha=layer.GetAlphaRef(ci);
 425        if (!fTrackToFollow->Propagate(alpha, x, GetBz())) {
 426          //Warning("TakeNextProlongation","propagation failed !\n");
 427           continue;
 428        }
 429        dz=7*TMath::Sqrt(fTrackToFollow->GetSigmaZ2() + kSigma2);
 430        dy=7*TMath::Sqrt(fTrackToFollow->GetSigmaY2() + kSigma2);
 431        z=fTrackToFollow->GetZ();
 432        y=fTrackToFollow->GetY();
 433     }
 434
 435     //if (TMath::Abs(fTrackToFollow.GetZ()-GetZ())>layer.GetR()+dz) continue;
 436
 437     if (TMath::Abs(z - c->GetZ()) > dz) continue;
 438     if (TMath::Abs(y - c->GetY()) > dy) continue;
 439
 440     Double_t ch2=fTrackToFollow->GetPredictedChi2(c);
 441     if (ch2 > kmaxChi2PerCluster) continue;
 442     chi2=ch2;
 443     cc=c; cci=ci;
 444     break;
 445   }
 446
 447   if (!cc) return 0;
 448
 449   if (!fTrackToFollow->Update(cc,chi2,(fI<<28)+cci)) {
 450      //Warning("TakeNextProlongation","filtering failed !\n");
 451      return 0;
 452   }
 453   Double_t xx0 = (fI > 2) ? 0.008 : 0.003;  // Rough layer thickness
 454   Double_t x0  = 9.36; // Radiation length of Si [cm]
 455   Double_t rho = 2.33; // Density of Si [g/cm^3]
 456   Double_t mass = fTrackToFollow->GetMass();
 457   fTrackToFollow->CorrectForMeanMaterial(xx0, xx0*x0*rho, mass, kTRUE);
 458   layer.ResetTrack(*fTrackToFollow);
 459
 460
 461   return 1;
 462 }
 463
 464 Int_t AliITSUTrackerCooked::PropagateBack(AliESDEvent *event) {
 465   //--------------------------------------------------------------------
 466   // Here, we implement the Kalman smoother ?
 467   // The clusters must already be loaded
 468   //--------------------------------------------------------------------
 469   Int_t n=event->GetNumberOfTracks();
 470   Int_t ntrk=0;
 471   Int_t ngood=0;
 472   for (Int_t i=0; i<n; i++) {
 473       AliESDtrack *esdTrack=event->GetTrack(i);
 474
 475       if (!esdTrack->IsOn(AliESDtrack::kITSin)) continue;
 476       if ( esdTrack->IsOn(AliESDtrack::kTPCin)) continue;//skip a TPC+ITS track
 477
 478       AliITSUTrackCooked track(*esdTrack);
 479
 480       ResetTrackToFollow(track);
 481
 482       fTrackToFollow->ResetCovariance(10.); fTrackToFollow->ResetClusters();
 483       if (RefitAt(40., fTrackToFollow, &track)) {
 484
 485          CookLabel(fTrackToFollow, 0.); //For comparison only
 486          Int_t label=fTrackToFollow->GetLabel();
 487          if (label>0) ngood++;
 488
 489          esdTrack->UpdateTrackParams(fTrackToFollow,AliESDtrack::kITSout);
 490          //UseClusters(fTrackToFollow);
 491          ntrk++;
 492       }
 493   }
 494
 495   Info("PropagateBack","Back propagated tracks: %d",ntrk);
 496   if (ntrk)
 497   Info("PropagateBack","Good tracks/back propagated: %f",Float_t(ngood)/ntrk);
 498
 499   if (!fSAonly) AliITSUTrackerGlo::PropagateBack(event);
 500
 501   return 0;
 502 }
 503
 504 Bool_t AliITSUTrackerCooked::
 505 RefitAt(Double_t xx, AliITSUTrackCooked *t, const AliITSUTrackCooked *c) {
 506   //--------------------------------------------------------------------
 507   // This function refits the track "t" at the position "x" using
 508   // the clusters from "c"
 509   //--------------------------------------------------------------------
 510   Int_t index[kNLayers];
 511   Int_t k;
 512   for (k=0; k<kNLayers; k++) index[k]=-1;
 513   Int_t nc=c->GetNumberOfClusters();
 514   for (k=0; k<nc; k++) {
 515     Int_t idx=c->GetClusterIndex(k), nl=(idx&0xf0000000)>>28;
 516     index[nl]=idx;
 517   }
 518
 519   Int_t from, to, step;
 520   if (xx > t->GetX()) {
 521       from=0; to=kNLayers;
 522       step=+1;
 523   } else {
 524       from=kNLayers-1; to=-1;
 525       step=-1;
 526   }
 527
 528   for (Int_t i=from; i != to; i += step) {
 529      Int_t idx=index[i];
 530      if (idx>=0) {
 531         const AliCluster *cl=GetCluster(idx);
 532         Float_t xr,ar; cl->GetXAlphaRefPlane(xr, ar);
 533         if (!t->Propagate(Double_t(ar), Double_t(xr), GetBz())) {
 534            //Warning("RefitAt","propagation failed !\n");
 535            return kFALSE;
 536         }
 537         Double_t chi2=t->GetPredictedChi2(cl);
 538         if (chi2 < kmaxChi2PerCluster) t->Update(cl, chi2, idx);
 539      } else {
 540         Double_t r=fgLayers[i].GetR();
 541         Double_t phi,z;
 542         if (!t->GetPhiZat(r,phi,z)) {
 543            //Warning("RefitAt","failed to estimate track !\n");
 544            return kFALSE;
 545         }
 546         if (!t->Propagate(phi, r, GetBz())) {
 547            //Warning("RefitAt","propagation failed !\n");
 548            return kFALSE;
 549         }
 550      }
 551      Double_t xx0 = (i > 2) ? 0.008 : 0.003;  // Rough layer thickness
 552      Double_t x0  = 9.36; // Radiation length of Si [cm]
 553      Double_t rho = 2.33; // Density of Si [g/cm^3]
 554      Double_t mass = t->GetMass();
 555      t->CorrectForMeanMaterial(xx0, -step*xx0*x0*rho, mass, kTRUE);
 556   }
 557
 558   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
 559   return kTRUE;
 560 }
 561
 562 Int_t AliITSUTrackerCooked::RefitInward(AliESDEvent *event) {
 563   //--------------------------------------------------------------------
 564   // Some final refit, after the outliers get removed by the smoother ?
 565   // The clusters must be loaded
 566   //--------------------------------------------------------------------
 567   Int_t n=event->GetNumberOfTracks();
 568   Int_t ntrk=0;
 569   Int_t ngood=0;
 570   for (Int_t i=0; i<n; i++) {
 571       AliESDtrack *esdTrack=event->GetTrack(i);
 572
 573       if (!esdTrack->IsOn(AliESDtrack::kITSout)) continue;
 574       if ( esdTrack->IsOn(AliESDtrack::kTPCin)) continue;//skip a TPC+ITS track
 575
 576       AliITSUTrackCooked track(*esdTrack);
 577       ResetTrackToFollow(track);
 578
 579       fTrackToFollow->ResetCovariance(10.); fTrackToFollow->ResetClusters();
 580       if (!RefitAt(2.1, fTrackToFollow, &track)) continue;
 581       //Cross the beam pipe
 582       if (!fTrackToFollow->PropagateTo(1.8, 2.27e-3, 35.28*1.848)) continue;
 583
 584       CookLabel(fTrackToFollow, 0.); //For comparison only
 585       Int_t label=fTrackToFollow->GetLabel();
 586       if (label>0) ngood++;
 587
 588       esdTrack->UpdateTrackParams(fTrackToFollow,AliESDtrack::kITSrefit);
 589       //esdTrack->RelateToVertex(event->GetVertex(),GetBz(),33.);
 590       //UseClusters(fTrackToFollow);
 591       ntrk++;
 592   }
 593
 594   Info("RefitInward","Refitted tracks: %d",ntrk);
 595   if (ntrk)
 596   Info("RefitInward","Good tracks/refitted: %f",Float_t(ngood)/ntrk);
 597
 598   if (!fSAonly) AliITSUTrackerGlo::RefitInward(event);
 599
 600   return 0;
 601 }
 602
 603 Int_t AliITSUTrackerCooked::LoadClusters(TTree *cTree) {
 604   //--------------------------------------------------------------------
 605   // This function reads the ITSU clusters from the tree,
 606   // sort them, distribute over the internal tracker arrays, etc
 607   //--------------------------------------------------------------------
 608   if (!cTree) {
 609      AliFatal("No cluster tree !");
 610      return 1;
 611   }
 612
 613   AliITSUTrackerGlo::LoadClusters(cTree);
 614
 615   for (Int_t i=0; i<kNLayers; i++) {
 616       TClonesArray *clusters=fReconstructor->GetClusters(i);
 617       switch (i) {
 618       case kSeedingLayer1:
 619       case kSeedingLayer2:
 620       case kSeedingLayer3:
 621          fgLayers[i].InsertClusters(clusters,kTRUE,fSAonly);
 622          break;
 623       default:
 624          fgLayers[i].InsertClusters(clusters,kFALSE,fSAonly);
 625          break;
 626       }
 627   }
 628
 629   return 0;
 630 }
 631
 632 void AliITSUTrackerCooked::UnloadClusters() {
 633   //--------------------------------------------------------------------
 634   // This function unloads ITSU clusters from the RAM
 635   //--------------------------------------------------------------------
 636   AliITSUTrackerGlo::UnloadClusters();
 637   for (Int_t i=0; i<kNLayers; i++) fgLayers[i].DeleteClusters();
 638 }
 639
 640 AliCluster *AliITSUTrackerCooked::GetCluster(Int_t index) const {
 641   //--------------------------------------------------------------------
 642   //       Return pointer to a given cluster
 643   //--------------------------------------------------------------------
 644     Int_t l=(index & 0xf0000000) >> 28;
 645     Int_t c=(index & 0x0fffffff) >> 00;
 646     return fgLayers[l].GetCluster(c);
 647 }
 648
 649 AliITSUTrackerCooked::AliITSUlayer::AliITSUlayer():
 650   fR(0),
 651   fN(0),
 652   fNsel(0),
 653   fI(0),
 654   fTrack(0)
 655 {
 656   //--------------------------------------------------------------------
 657   // This default constructor needs to be provided
 658   //--------------------------------------------------------------------
 659   for (Int_t i=0; i<kMaxClusterPerLayer; i++) fClusters[i]=0;
 660   for (Int_t i=0; i<kMaxSelected; i++) fIndex[i]=-1;
 661 }
 662
 663 AliITSUTrackerCooked::AliITSUlayer::~AliITSUlayer()
 664 {
 665   //--------------------------------------------------------------------
 666   // Simple destructor
 667   //--------------------------------------------------------------------
 668   delete fTrack;
 669 }
 670
 671 void
 672 AliITSUTrackerCooked::AliITSUlayer::ResetTrack(const AliITSUTrackCooked &t) {
 673   //--------------------------------------------------------------------
 674   // Replace the track estimate at this layer
 675   //--------------------------------------------------------------------
 676    delete fTrack;
 677    fTrack=new AliITSUTrackCooked(t);
 678 }
 679
 680 void AliITSUTrackerCooked::AliITSUlayer::
 681 InsertClusters(TClonesArray *clusters, Bool_t seedingLayer, Bool_t saOnly)
 682 {
 683   //--------------------------------------------------------------------
 684   // Load clusters to this layer
 685   //--------------------------------------------------------------------
 686   Int_t ncl=clusters->GetEntriesFast();
 687   Double_t r=0.;
 688   for (Int_t i=0; i<ncl; i++) {
 689      AliITSUClusterPix *c=(AliITSUClusterPix*)clusters->UncheckedAt(i);
 690      if (!saOnly) if (c->IsClusterUsed()) continue;
 691      c->GoToFrameGlo();
 692      Double_t x=c->GetX(), y=c->GetY();
 693      r += TMath::Sqrt(x*x + y*y);
 694      if (!seedingLayer) c->GoToFrameTrk();
 695      //if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 696      if (InsertCluster(c)) break;
 697   }
 698   if (fN) fR = r/fN;
 699   const Float_t pi2 = 2.*TMath::Pi();
 700   for (Int_t i=0; i<fN; i++) {
 701       AliCluster *c=fClusters[i];
 702       c->GetXAlphaRefPlane(fXRef[i],fAlphaRef[i]);
 703       Float_t xyz[3]; c->GetGlobalXYZ(xyz);
 704       Float_t phi=TMath::ATan2(xyz[1],xyz[0]);
 705       if (phi<0.) phi+=pi2;
 706       else if (phi >= pi2) phi-=pi2;
 707       fPhi[i]=phi;
 708   }
 709 }
 710
 711 void AliITSUTrackerCooked::AliITSUlayer::DeleteClusters()
 712 {
 713   //--------------------------------------------------------------------
 714   // Load clusters to this layer
 715   //--------------------------------------------------------------------
 716   //for (Int_t i=0; i<fN; i++) {delete fClusters[i]; fClusters[i]=0;}
 717   fN=0;
 718 }
 719
 720 Int_t
 721 AliITSUTrackerCooked::AliITSUlayer::InsertCluster(AliCluster *c) {
 722   //--------------------------------------------------------------------
 723   // This function inserts a cluster to this layer in increasing
 724   // order of the cluster's fZ
 725   //--------------------------------------------------------------------
 726   if (fN>=kMaxClusterPerLayer) {
 727      ::Error("InsertCluster","Too many clusters !\n");
 728      return 1;
 729   }
 730   if (fN==0) fClusters[0]=c;
 731   else {
 732      Int_t i=FindClusterIndex(c->GetZ());
 733      Int_t k=fN-i;
 734      memmove(fClusters+i+1 ,fClusters+i,k*sizeof(AliCluster*));
 735      fClusters[i]=c;
 736   }
 737   fN++;
 738   return 0;
 739 }
 740
 741 Int_t
 742 AliITSUTrackerCooked::AliITSUlayer::FindClusterIndex(Double_t z) const {
 743   //--------------------------------------------------------------------
 744   // This function returns the index of the first
 745   // with its fZ >= "z".
 746   //--------------------------------------------------------------------
 747   if (fN==0) return 0;
 748
 749   Int_t b=0;
 750   if (z <= fClusters[b]->GetZ()) return b;
 751
 752   Int_t e=b+fN-1;
 753   if (z > fClusters[e]->GetZ()) return e+1;
 754
 755   Int_t m=(b+e)/2;
 756   for (; b<e; m=(b+e)/2) {
 757     if (z > fClusters[m]->GetZ()) b=m+1;
 758     else e=m;
 759   }
 760   return m;
 761 }
 762
 763 void AliITSUTrackerCooked::AliITSUlayer::
 764 SelectClusters(Float_t phi, Float_t dy, Float_t z, Float_t dz) {
 765   //--------------------------------------------------------------------
 766   // This function selects clusters within the "road"
 767   //--------------------------------------------------------------------
 768   fNsel=0;
 769
 770   Float_t dphi=dy/fR;
 771   Float_t phiMin=phi-dphi;
 772   Float_t phiMax=phi+dphi;
 773   Float_t zMin=z-dz;
 774   Float_t zMax=z+dz;
 775
 776   Int_t i=FindClusterIndex(zMin), imax=FindClusterIndex(zMax);
 777   for (; i<imax; i++) {
 778       Float_t cphi=fPhi[i];
 779       if (cphi <= phiMin) continue;
 780       if (cphi >  phiMax) continue;
 781       AliCluster *c=fClusters[i];
 782       if (c->IsClusterUsed()) continue;
 783
 784       fIndex[fNsel++]=i;
 785       if (fNsel>=kMaxSelected) break;
 786   }
 787 }
 788
 789 const AliCluster *AliITSUTrackerCooked::AliITSUlayer::GetNextCluster(Int_t &ci){
 790   //--------------------------------------------------------------------
 791   // This function returns clusters within the "road"
 792   //--------------------------------------------------------------------
 793   if (fI<fNsel) {
 794       ci=fIndex[fI++];
 795       return fClusters[ci];
 796   }
 797   ci=-1;
 798   return 0;
 799 }
 800