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