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