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 <Riostream.h>
   7 #include <TTree.h>
   8 #include <TClonesArray.h>
   9
  10 #include "AliLog.h"
  11 #include "AliESDEvent.h"
  12 #include "AliITSUClusterPix.h"
  13 #include "AliITSUGeomTGeo.h"
  14 #include "AliITSUTrackerCooked.h"
  15 #include "AliITSUTrackCooked.h"
  16
  17 using std::cout;
  18 using std::endl;
  19
  20 ClassImp(AliITSUTrackerCooked)
  21
  22 //************************************************
  23 // Constants hardcoded for the moment:
  24 //************************************************
  25 // radial positions of layers: default contructor
  26 const
  27 Double_t klRadius[7]={2.34, 3.15, 3.93, 19.61, 24.55, 34.39, 39.34}; //tdr6
  28 // seed "windows" in z and phi: MakeSeeds
  29 const Double_t kzWin=0.33, kpWin=3.14/4;
  30 // Maximal accepted impact parameters for the seeds
  31 const Double_t kmaxDCAxy=3.;
  32 const Double_t kmaxDCAz= 3.;
  33 // Layers for the seeding
  34 const Int_t kSeedingLayer1=6, kSeedingLayer2=4, kSeedingLayer3=5;
  35
  36
  37 //************************************************
  38 // TODO:
  39 //************************************************
  40 // Seeding:
  41 // Precalculate cylidnrical (r,phi) for the clusters;
  42 // use exact r's for the clusters
  43
  44
  45 AliITSUTrackerCooked::AliITSUlayer
  46               AliITSUTrackerCooked::fgLayers[AliITSUTrackerCooked::kNLayers];
  47
  48 AliITSUTrackerCooked::AliITSUlayer::AliITSUlayer():
  49   fR(0),
  50   fN(0),
  51   fNsel(0)
  52 {
  53   //--------------------------------------------------------------------
  54   // This default constructor needs to be provided
  55   //--------------------------------------------------------------------
  56   for (Int_t i=0; i<kMaxClusterPerLayer; i++) {
  57     fClusters[i]=0;
  58     fIndex[i]=-1;
  59   }
  60 }
  61
  62 void AliITSUTrackerCooked::
  63   AliITSUlayer::InsertClusters(TClonesArray *clusters, Bool_t seedingLayer)
  64 {
  65   //--------------------------------------------------------------------
  66   // Load clusters to this layer
  67   //--------------------------------------------------------------------
  68   Int_t ncl=clusters->GetEntriesFast();
  69
  70   while (ncl--) {
  71      AliITSUClusterPix *c=(AliITSUClusterPix*)clusters->UncheckedAt(ncl);
  72      (seedingLayer) ? c->GoToFrameGlo() : c->GoToFrameTrk();
  73      //if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
  74      InsertCluster(new AliITSUClusterPix(*c));
  75   }
  76
  77 }
  78
  79 void AliITSUTrackerCooked::AliITSUlayer::DeleteClusters()
  80 {
  81   //--------------------------------------------------------------------
  82   // Load clusters to this layer
  83   //--------------------------------------------------------------------
  84   for (Int_t i=0; i<fN; i++) delete fClusters[i];
  85   fN=0;
  86 }
  87
  88 Int_t
  89 AliITSUTrackerCooked::AliITSUlayer::InsertCluster(AliCluster *c) {
  90   //--------------------------------------------------------------------
  91   // This function inserts a cluster to this layer in increasing
  92   // order of the cluster's fZ
  93   //--------------------------------------------------------------------
  94   if (fN>=kMaxClusterPerLayer) {
  95      ::Error("InsertCluster","Too many clusters !\n");
  96      return 1;
  97   }
  98   if (fN==0) fClusters[0]=c;
  99   else {
 100      Int_t i=FindClusterIndex(c->GetZ());
 101      Int_t k=fN-i;
 102      memmove(fClusters+i+1 ,fClusters+i,k*sizeof(AliCluster*));
 103      fClusters[i]=c;
 104   }
 105   fN++;
 106   return 0;
 107 }
 108
 109 Int_t
 110 AliITSUTrackerCooked::AliITSUlayer::FindClusterIndex(Double_t z) const {
 111   //--------------------------------------------------------------------
 112   // This function returns the index of the first
 113   // with its fZ >= "z".
 114   //--------------------------------------------------------------------
 115   if (fN==0) return 0;
 116
 117   Int_t b=0;
 118   if (z <= fClusters[b]->GetZ()) return b;
 119
 120   Int_t e=b+fN-1;
 121   if (z > fClusters[e]->GetZ()) return e+1;
 122
 123   Int_t m=(b+e)/2;
 124   for (; b<e; m=(b+e)/2) {
 125     if (z > fClusters[m]->GetZ()) b=m+1;
 126     else e=m;
 127   }
 128   return m;
 129 }
 130
 131
 132 AliITSUTrackerCooked::AliITSUTrackerCooked():
 133   AliTracker(),fSeeds(0)
 134 {
 135   //--------------------------------------------------------------------
 136   // This default constructor needs to be provided
 137   //--------------------------------------------------------------------
 138
 139   AliITSUGeomTGeo *gm  = new AliITSUGeomTGeo(kTRUE,kTRUE);
 140   AliITSUClusterPix::SetGeom(gm);
 141
 142   for (Int_t i=0; i<kNLayers; i++) fgLayers[i].SetR(klRadius[i]);
 143
 144   // Some default primary vertex
 145   Double_t xyz[]={0.,0.,0.};
 146   Double_t ers[]={2.,2.,2.};
 147
 148   SetVertex(xyz,ers);
 149
 150 }
 151
 152 AliITSUTrackerCooked::~AliITSUTrackerCooked()
 153 {
 154   //--------------------------------------------------------------------
 155   // Virtual destructor
 156   //--------------------------------------------------------------------
 157
 158   UnloadClusters();
 159
 160   if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
 161   fSeeds=0;
 162 }
 163
 164 static Double_t
 165 f1(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
 166 {
 167     //-----------------------------------------------------------------
 168     // Initial approximation of the track curvature
 169     //-----------------------------------------------------------------
 170     Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 171     Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 172                     (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 173     Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 174                     (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 175
 176     Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 177
 178     return -xr*yr/sqrt(xr*xr+yr*yr);
 179 }
 180
 181 static Double_t
 182 f2(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3)
 183 {
 184     //-----------------------------------------------------------------
 185     // Initial approximation of the track curvature times center of curvature
 186     //-----------------------------------------------------------------
 187     Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
 188     Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
 189                     (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
 190     Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
 191                     (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
 192
 193     Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
 194
 195     return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
 196 }
 197
 198 static Double_t
 199 f3(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t z1, Double_t z2)
 200 {
 201     //-----------------------------------------------------------------
 202     // Initial approximation of the tangent of the track dip angle
 203     //-----------------------------------------------------------------
 204     return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
 205 }
 206
 207 Bool_t AliITSUTrackerCooked::
 208 AddCookedSeed(const Float_t r1[3], Int_t l1, Int_t i1,
 209               const Float_t r2[3], Int_t l2, Int_t i2,
 210               const AliCluster *c3,Int_t l3, Int_t i3)
 211 {
 212     //--------------------------------------------------------------------
 213     // This is the main cooking function.
 214     // Creates seed parameters out of provided clusters.
 215     //--------------------------------------------------------------------
 216     Float_t x,a;
 217     if (!c3->GetXAlphaRefPlane(x,a)) return kFALSE;
 218
 219     Double_t ca=TMath::Cos(a), sa=TMath::Sin(a);
 220     Double_t x1 = r1[0]*ca + r1[1]*sa,
 221              y1 =-r1[0]*sa + r1[1]*ca, z1 = r1[2];
 222     Double_t x2 = r2[0]*ca + r2[1]*sa,
 223              y2 =-r2[0]*sa + r2[1]*ca, z2 = r2[2];
 224     Double_t x3 = x,  y3 = c3->GetY(), z3 = c3->GetZ();
 225
 226     Double_t par[5];
 227     par[0]=y3;
 228     par[1]=z3;
 229     Double_t crv=f1(x1, y1, x2, y2, x3, y3); //curvature
 230     Double_t cx0=f2(x1, y1, x2, y2, x3, y3); //curvature*x0
 231     Double_t tgl12=f3(x1, y1, x2, y2, z1, z2);
 232     Double_t tgl23=f3(x2, y2, x3, y3, z2, z3);
 233
 234     Double_t sf=x*crv - cx0;
 235     if (TMath::Abs(sf) >= kAlmost1) return kFALSE;
 236     par[2]=sf;
 237
 238     par[3]=0.5*(tgl12 + tgl23);
 239     Double_t bz=GetBz();
 240     par[4]=(TMath::Abs(bz) < kAlmost0Field) ? kAlmost0 : crv/(bz*kB2C);
 241
 242     Double_t cov[15];
 243     for (Int_t i=0; i<15; i++) cov[i]=0.;
 244     cov[0] =0.0005*0.0005*10;
 245     cov[2] =0.0005*0.0005*10;
 246     cov[5] =0.007*0.007*10;   //FIXME all these lines
 247     cov[9] =0.007*0.007*10;
 248     cov[14]=0.1*0.1*10;
 249
 250     AliITSUTrackCooked *seed=new AliITSUTrackCooked();
 251     seed->Set(Double_t(x), Double_t(a), par, cov);
 252
 253     Float_t dz[2];
 254     seed->GetDZ(GetX(),GetY(),GetZ(),GetBz(),dz);
 255     if (TMath::Abs(dz[0]) > kmaxDCAxy) {delete seed; return kFALSE;}
 256     if (TMath::Abs(dz[1]) > kmaxDCAz ) {delete seed; return kFALSE;}
 257
 258     seed->SetClusterIndex(l1,i1);
 259     seed->SetClusterIndex(l2,i2);
 260     seed->SetClusterIndex(l3,i3);
 261
 262     fSeeds->AddLast(seed);
 263
 264     return kTRUE;
 265 }
 266
 267 Int_t AliITSUTrackerCooked::MakeSeeds() {
 268   //--------------------------------------------------------------------
 269   // This is the main pattern recongition function.
 270   // Creates seeds out of two clusters and another point.
 271   //--------------------------------------------------------------------
 272    if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
 273    fSeeds=new TObjArray(77777);
 274
 275    Double_t zv=GetZ();
 276
 277    AliITSUlayer &layer1=fgLayers[kSeedingLayer1];
 278    AliITSUlayer &layer2=fgLayers[kSeedingLayer2];
 279    AliITSUlayer &layer3=fgLayers[kSeedingLayer3];
 280    Double_t r1=layer1.GetR();
 281    Double_t r2=layer2.GetR();
 282    Double_t r3=layer3.GetR();
 283
 284    Int_t nClusters1=layer1.GetNumberOfClusters();
 285    Int_t nClusters2=layer2.GetNumberOfClusters();
 286    Int_t nClusters3=layer3.GetNumberOfClusters();
 287    for (Int_t n1=0; n1<nClusters1; n1++) {
 288      AliCluster *c1=layer1.GetCluster(n1);
 289      //
 290      //Int_t lab=c1->GetLabel(0);
 291      //
 292      Double_t z1=c1->GetZ();
 293      Float_t xyz1[3]; c1->GetGlobalXYZ(xyz1);
 294      Double_t phi1=TMath::ATan2(xyz1[1],xyz1[0]);
 295      Double_t zr2=zv + r2/r1*(z1-zv);
 296      Int_t start2=layer2.FindClusterIndex(zr2-kzWin);
 297      for (Int_t n2=start2; n2<nClusters2; n2++) {
 298          AliCluster *c2=layer2.GetCluster(n2);
 299          //
 300          //if (c2->GetLabel(0)!=lab) continue;
 301          //
 302          Double_t z2=c2->GetZ();
 303
 304          if (z2 > (zr2+kzWin)) break;  //check in Z
 305          Float_t xyz2[3]; c2->GetGlobalXYZ(xyz2);
 306          Double_t phi2=TMath::ATan2(xyz2[1],xyz2[0]);
 307          if (TMath::Abs(phi2-phi1) > kpWin) continue;  //check in Phi
 308
 309          Double_t zr3=z1 + (r3-r1)/(r2-r1)*(z2-z1);
 310          Double_t d13=r1-r3, d32=r3-r2, d=r1-r2;
 311          Double_t phir3=d32/d*phi1 + d13/d*phi2;  // FIXME
 312          Int_t start3=layer3.FindClusterIndex(zr3-kzWin/2);
 313          for (Int_t n3=start3; n3<nClusters3; n3++) {
 314              AliCluster *c3=layer3.GetCluster(n3);
 315              //
 316              //if (c3->GetLabel(0)!=lab) continue;
 317              //
 318              Double_t z3=c3->GetZ();
 319
 320              if (z3 > (zr3+kzWin/2)) break;  //check in Z
 321              Float_t xyz3[3]; c3->GetGlobalXYZ(xyz3);
 322              Double_t phi3=TMath::ATan2(xyz3[1],xyz3[0]);
 323              if (TMath::Abs(phir3-phi3) > kpWin/100) continue;  //check in Phi
 324
 325              AliITSUClusterPix cc(*((AliITSUClusterPix*)c2));
 326              cc.GoToFrameTrk();
 327              AddCookedSeed(xyz1, kSeedingLayer1, n1,
 328                            xyz3, kSeedingLayer3, n3,
 329                            &cc,  kSeedingLayer2, n2);
 330
 331          }
 332      }
 333    }
 334    fSeeds->Sort();
 335    return fSeeds->GetEntriesFast();
 336 }
 337
 338 Int_t AliITSUTrackerCooked::Clusters2Tracks(AliESDEvent *event) {
 339   //--------------------------------------------------------------------
 340   // This is the main tracking function
 341   // The clusters must already be loaded
 342   //--------------------------------------------------------------------
 343
 344   Int_t nSeeds=MakeSeeds();
 345
 346   // Possibly, icrement the seeds with additional clusters (Kalman)
 347
 348   // Possibly, (re)fit the found tracks
 349
 350
 351   Int_t nClusters1=fgLayers[kSeedingLayer1].GetNumberOfClusters();
 352   Int_t nClusters2=fgLayers[kSeedingLayer2].GetNumberOfClusters();
 353   cout<<nClusters1<<' '<<nClusters2<<endl;
 354
 355     Int_t ngood=0;
 356     for (Int_t s=0; s<nSeeds; s++) {
 357         AliITSUTrackCooked *track=(AliITSUTrackCooked*)fSeeds->At(s);
 358         CookLabel(track,0.);
 359         if (track->GetLabel()>0) {ngood++;
 360             //cout<<track->GetLabel()<<endl;
 361         }
 362         AliESDtrack iotrack;
 363         iotrack.UpdateTrackParams(track,AliESDtrack::kITSin);
 364         iotrack.SetLabel(track->GetLabel());
 365         event->AddTrack(&iotrack);
 366     }
 367     cout<<"Good tracks "<<ngood<<endl;
 368     cout<<"Good tracks/seeds "<<Float_t(ngood)/nSeeds<<endl;
 369
 370   if (fSeeds) {fSeeds->Delete(); delete fSeeds;}
 371   fSeeds=0;
 372
 373   return 0;
 374 }
 375
 376 Int_t AliITSUTrackerCooked::PropagateBack(AliESDEvent *event) {
 377   //--------------------------------------------------------------------
 378   // Here, we implement the Kalman smoother ?
 379   // The clusters must already be loaded
 380   //--------------------------------------------------------------------
 381     Int_t n=event->GetNumberOfTracks();
 382     for (Int_t i=0; i<n; i++) {
 383         AliESDtrack *esdTrack=event->GetTrack(i);
 384         AliITSUTrackCooked track(*esdTrack);
 385         esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSout);
 386     }
 387
 388   return 0;
 389 }
 390
 391 Int_t AliITSUTrackerCooked::RefitInward(AliESDEvent *event) {
 392   //--------------------------------------------------------------------
 393   // Some final refit, after the outliers get removed by the smoother ?
 394   // The clusters must be loaded
 395   //--------------------------------------------------------------------
 396     Int_t n=event->GetNumberOfTracks();
 397     for (Int_t i=0; i<n; i++) {
 398         AliESDtrack *esdTrack=event->GetTrack(i);
 399         AliITSUTrackCooked track(*esdTrack);
 400         esdTrack->UpdateTrackParams(&track,AliESDtrack::kITSrefit);
 401     }
 402
 403   return 0;
 404 }
 405
 406 Int_t AliITSUTrackerCooked::LoadClusters(TTree *cTree) {
 407   //--------------------------------------------------------------------
 408   // This function reads the ITSU clusters from the tree,
 409   // sort them, distribute over the internal tracker arrays, etc
 410   //--------------------------------------------------------------------
 411   if (!cTree) {
 412      AliFatal("No cluster tree !");
 413      return 1;
 414   }
 415
 416   //This TClonesArray is not the owner of the clusters
 417   TClonesArray dummy("AliITSUClusterPix",kMaxClusterPerLayer), *clusters=&dummy;
 418
 419   for (Int_t i=0; i<kNLayers; i++) {
 420       TBranch *br = cTree->GetBranch(Form("ITSRecPoints%d",i));
 421       if (!br) AliFatal(Form("No cluster branch for layer %d",i));
 422       br->SetAddress(&clusters);
 423       br->GetEvent(0);
 424       switch (i) {
 425       case kSeedingLayer1:
 426       case kSeedingLayer2:
 427       case kSeedingLayer3:
 428          fgLayers[i].InsertClusters(clusters,kTRUE);
 429          break;
 430       default:
 431          fgLayers[i].InsertClusters(clusters,kFALSE);
 432          break;
 433       }
 434       clusters->Delete();
 435   }
 436
 437   return 0;
 438 }
 439
 440 void AliITSUTrackerCooked::UnloadClusters() {
 441   //--------------------------------------------------------------------
 442   // This function unloads ITSU clusters from the RAM
 443   //--------------------------------------------------------------------
 444   for (Int_t i=0; i<kNLayers; i++) fgLayers[i].DeleteClusters();
 445 }
 446
 447 AliCluster *AliITSUTrackerCooked::GetCluster(Int_t index) const {
 448   //--------------------------------------------------------------------
 449   //       Return pointer to a given cluster
 450   //--------------------------------------------------------------------
 451     Int_t l=(index & 0xf0000000) >> 28;
 452     Int_t c=(index & 0x0fffffff) >> 00;
 453     return fgLayers[l].GetCluster(c);
 454 }