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