ITS/AliITStrackerV2.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-------------------------------------------------------------------------
  17 //               Implementation of the ITS tracker class
  18 //    It reads AliITSRecPoint clusters and creates AliITStrackV2 tracks
  19 //                   and fills with them the ESD
  20 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
  21 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  22 //-------------------------------------------------------------------------
  23
  24 #include <new>
  25
  26 #include <TFile.h>
  27 #include <TTree.h>
  28 #include <TRandom.h>
  29
  30 #include "AliITSgeom.h"
  31 #include "AliITSRecPoint.h"
  32 #include "AliESD.h"
  33 #include "AliITSRecPoint.h"
  34 #include "AliITStrackerV2.h"
  35
  36 ClassImp(AliITStrackerV2)
  37
  38 AliITStrackerV2::AliITSlayer AliITStrackerV2::fgLayers[kMaxLayer]; // ITS layers
  39
  40 AliITStrackerV2::AliITStrackerV2(const AliITSgeom *geom) : AliTracker() {
  41   //--------------------------------------------------------------------
  42   //This is the AliITStrackerV2 constructor
  43   //--------------------------------------------------------------------
  44   AliITSgeom *g=(AliITSgeom*)geom;
  45
  46   Float_t x,y,z;  Int_t i;
  47   for (i=1; i<kMaxLayer+1; i++) {
  48     Int_t nlad=g->GetNladders(i);
  49     Int_t ndet=g->GetNdetectors(i);
  50
  51     g->GetTrans(i,1,1,x,y,z);
  52     Double_t r=TMath::Sqrt(x*x + y*y);
  53     Double_t poff=TMath::ATan2(y,x);
  54     Double_t zoff=z;
  55
  56     g->GetTrans(i,1,2,x,y,z);
  57     r += TMath::Sqrt(x*x + y*y);
  58     g->GetTrans(i,2,1,x,y,z);
  59     r += TMath::Sqrt(x*x + y*y);
  60     g->GetTrans(i,2,2,x,y,z);
  61     r += TMath::Sqrt(x*x + y*y);
  62     r*=0.25;
  63
  64     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
  65
  66     for (Int_t j=1; j<nlad+1; j++) {
  67       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
  68         Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift);
  69         Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
  70
  71         Double_t phi=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
  72         phi+=TMath::Pi()/2;
  73         if (i==1) phi+=TMath::Pi();
  74
  75         if (phi<0) phi+=TMath::TwoPi();
  76         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
  77
  78         Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
  79         Double_t r=x*cp+y*sp;
  80
  81         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
  82         new(&det) AliITSdetector(r,phi);
  83       }
  84     }
  85
  86   }
  87
  88   fI=kMaxLayer;
  89
  90   fPass=0;
  91   fConstraint[0]=1; fConstraint[1]=0;
  92
  93   Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV};
  94   SetVertex(xyz,ers);
  95
  96   for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=kLayersNotToSkip[i];
  97   fLastLayerToTrackTo=kLastLayerToTrackTo;
  98
  99 }
 100
 101 void AliITStrackerV2::SetLayersNotToSkip(Int_t *l) {
 102   //--------------------------------------------------------------------
 103   //This function set masks of the layers which must be not skipped
 104   //--------------------------------------------------------------------
 105   for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=l[i];
 106 }
 107
 108 Int_t AliITStrackerV2::LoadClusters(TTree *cTree) {
 109   //--------------------------------------------------------------------
 110   //This function loads ITS clusters
 111   //--------------------------------------------------------------------
 112   TBranch *branch=cTree->GetBranch("ITSRecPoints");
 113   if (!branch) {
 114     Error("LoadClusters"," can't get the branch !\n");
 115     return 1;
 116   }
 117
 118   TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 119   branch->SetAddress(&clusters);
 120
 121   Int_t j=0;
 122   for (Int_t i=0; i<kMaxLayer; i++) {
 123     Int_t ndet=fgLayers[i].GetNdetectors();
 124     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 125
 126     Double_t r=fgLayers[i].GetR();
 127     Double_t circ=TMath::TwoPi()*r;
 128
 129     for (; j<jmax; j++) {
 130       if (!cTree->GetEvent(j)) continue;
 131       Int_t ncl=clusters->GetEntriesFast();
 132       while (ncl--) {
 133         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 134
 135         Int_t idx=c->GetDetectorIndex();
 136         Double_t y=r*fgLayers[i].GetDetector(idx).GetPhi()+c->GetY();
 137         if (y>circ) y-=circ; else if (y<0) y+=circ;
 138         c->SetPhiR(y);
 139
 140         fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 141       }
 142       clusters->Delete();
 143     }
 144     fgLayers[i].ResetRoad(); //road defined by the cluster density
 145   }
 146
 147   return 0;
 148 }
 149
 150 void AliITStrackerV2::UnloadClusters() {
 151   //--------------------------------------------------------------------
 152   //This function unloads ITS clusters
 153   //--------------------------------------------------------------------
 154   for (Int_t i=0; i<kMaxLayer; i++) fgLayers[i].ResetClusters();
 155 }
 156
 157 static Int_t CorrectForDeadZoneMaterial(AliITStrackV2 *t) {
 158   //--------------------------------------------------------------------
 159   // Correction for the material between the TPC and the ITS
 160   // (should it belong to the TPC code ?)
 161   //--------------------------------------------------------------------
 162   Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ?
 163   Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
 164   Double_t yr=12.8, dr=0.03; // rods ?
 165   Double_t zm=0.2, dm=0.40;  // membrane
 166   //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
 167   Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)
 168
 169   if (t->GetX() > riw) {
 170      if (!t->PropagateTo(riw,diw,x0iw)) return 1;
 171      if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
 172      if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
 173      if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
 174      //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
 175      //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r);
 176      if (!t->PropagateTo(rs,ds)) return 1;
 177   } else if (t->GetX() < rs) {
 178      if (!t->PropagateTo(rs,-ds)) return 1;
 179      //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
 180      //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r);
 181      if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
 182      if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
 183   } else {
 184   ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
 185     return 1;
 186   }
 187
 188   return 0;
 189 }
 190
 191 Int_t AliITStrackerV2::Clusters2Tracks(AliESD *event) {
 192   //--------------------------------------------------------------------
 193   // This functions reconstructs ITS tracks
 194   // The clusters must be already loaded !
 195   //--------------------------------------------------------------------
 196   TObjArray itsTracks(15000);
 197
 198   {/* Read ESD tracks */
 199     Int_t nentr=event->GetNumberOfTracks();
 200     Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 201     while (nentr--) {
 202       AliESDtrack *esd=event->GetTrack(nentr);
 203
 204       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 205       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 206       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 207
 208       AliITStrackV2 *t=0;
 209       try {
 210         t=new AliITStrackV2(*esd);
 211       } catch (const Char_t *msg) {
 212         Warning("Clusters2Tracks",msg);
 213         delete t;
 214         continue;
 215       }
 216       if (TMath::Abs(t->GetD())>4) {
 217         delete t;
 218         continue;
 219       }
 220
 221       if (CorrectForDeadZoneMaterial(t)!=0) {
 222          Warning("Clusters2Tracks",
 223                  "failed to correct for the material in the dead zone !\n");
 224          delete t;
 225          continue;
 226       }
 227       itsTracks.AddLast(t);
 228     }
 229   } /* End Read ESD tracks */
 230
 231   itsTracks.Sort();
 232   Int_t nentr=itsTracks.GetEntriesFast();
 233
 234   Int_t ntrk=0;
 235   for (fPass=0; fPass<2; fPass++) {
 236      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 237      for (Int_t i=0; i<nentr; i++) {
 238        AliITStrackV2 *t=(AliITStrackV2*)itsTracks.UncheckedAt(i);
 239        if (t==0) continue;           //this track has been already tracked
 240        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 241
 242        ResetTrackToFollow(*t);
 243        ResetBestTrack();
 244
 245        for (FollowProlongation(); fI<kMaxLayer; fI++) {
 246           while (TakeNextProlongation()) FollowProlongation();
 247        }
 248
 249        if (fBestTrack.GetNumberOfClusters() == 0) continue;
 250
 251        if (fConstraint[fPass]) {
 252           ResetTrackToFollow(*t);
 253           if (!RefitAt(3.7, &fTrackToFollow, &fBestTrack)) continue;
 254           ResetBestTrack();
 255        }
 256
 257        fBestTrack.SetLabel(tpcLabel);
 258        fBestTrack.CookdEdx();
 259        CookLabel(&fBestTrack,0.); //For comparison only
 260        fBestTrack.UpdateESDtrack(AliESDtrack::kITSin);
 261        UseClusters(&fBestTrack);
 262        delete itsTracks.RemoveAt(i);
 263        ntrk++;
 264      }
 265   }
 266
 267   itsTracks.Delete();
 268
 269   Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
 270
 271   return 0;
 272 }
 273
 274 Int_t AliITStrackerV2::PropagateBack(AliESD *event) {
 275   //--------------------------------------------------------------------
 276   // This functions propagates reconstructed ITS tracks back
 277   // The clusters must be loaded !
 278   //--------------------------------------------------------------------
 279   Int_t nentr=event->GetNumberOfTracks();
 280   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 281
 282   Int_t ntrk=0;
 283   for (Int_t i=0; i<nentr; i++) {
 284      AliESDtrack *esd=event->GetTrack(i);
 285
 286      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 287      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 288
 289      AliITStrackV2 *t=0;
 290      try {
 291         t=new AliITStrackV2(*esd);
 292      } catch (const Char_t *msg) {
 293         Warning("PropagateBack",msg);
 294         delete t;
 295         continue;
 296      }
 297
 298      ResetTrackToFollow(*t);
 299
 300      // propagete to vertex [SR, GSI 17.02.2003]
 301      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
 302      if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {
 303        if (fTrackToFollow.PropagateToVertex()) {
 304           fTrackToFollow.StartTimeIntegral();
 305        }
 306        fTrackToFollow.PropagateTo(3.,-0.0028,65.19);
 307      }
 308
 309      fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
 310      if (RefitAt(49.,&fTrackToFollow,t)) {
 311         if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) {
 312           Warning("PropagateBack",
 313                   "failed to correct for the material in the dead zone !\n");
 314           delete t;
 315           continue;
 316         }
 317         fTrackToFollow.SetLabel(t->GetLabel());
 318         //fTrackToFollow.CookdEdx();
 319         CookLabel(&fTrackToFollow,0.); //For comparison only
 320         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 321         UseClusters(&fTrackToFollow);
 322         ntrk++;
 323      }
 324      delete t;
 325   }
 326
 327   Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
 328
 329   return 0;
 330 }
 331
 332 Int_t AliITStrackerV2::RefitInward(AliESD *event) {
 333   //--------------------------------------------------------------------
 334   // This functions refits ITS tracks using the
 335   // "inward propagated" TPC tracks
 336   // The clusters must be loaded !
 337   //--------------------------------------------------------------------
 338   Int_t nentr=event->GetNumberOfTracks();
 339   Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 340
 341   Int_t ntrk=0;
 342   for (Int_t i=0; i<nentr; i++) {
 343     AliESDtrack *esd=event->GetTrack(i);
 344
 345     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 346     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 347     if (esd->GetStatus()&AliESDtrack::kTPCout)
 348     if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 349
 350     AliITStrackV2 *t=0;
 351     try {
 352         t=new AliITStrackV2(*esd);
 353     } catch (const Char_t *msg) {
 354         Warning("RefitInward",msg);
 355         delete t;
 356         continue;
 357     }
 358
 359     if (CorrectForDeadZoneMaterial(t)!=0) {
 360        Warning("RefitInward",
 361                "failed to correct for the material in the dead zone !\n");
 362        delete t;
 363        continue;
 364     }
 365
 366     ResetTrackToFollow(*t);
 367     fTrackToFollow.ResetClusters();
 368
 369     //Refitting...
 370     if (RefitAt(3.7, &fTrackToFollow, t)) {
 371        fTrackToFollow.SetLabel(t->GetLabel());
 372        fTrackToFollow.CookdEdx();
 373        CookLabel(&fTrackToFollow,0.); //For comparison only
 374
 375        if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe
 376          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
 377          esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit);
 378          Float_t r[3]={0.,0.,0.};
 379          Double_t maxD=3.;
 380          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
 381          ntrk++;
 382        }
 383     }
 384     delete t;
 385   }
 386
 387   Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
 388
 389   return 0;
 390 }
 391
 392 AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
 393   //--------------------------------------------------------------------
 394   //       Return pointer to a given cluster
 395   //--------------------------------------------------------------------
 396   Int_t l=(index & 0xf0000000) >> 28;
 397   Int_t c=(index & 0x0fffffff) >> 00;
 398   return fgLayers[l].GetCluster(c);
 399 }
 400
 401
 402 void AliITStrackerV2::FollowProlongation() {
 403   //--------------------------------------------------------------------
 404   //This function finds a track prolongation
 405   //--------------------------------------------------------------------
 406   while (fI>fLastLayerToTrackTo) {
 407     Int_t i=fI-1;
 408
 409     AliITSlayer &layer=fgLayers[i];
 410     AliITStrackV2 &track=fTracks[i];
 411
 412     Double_t r=layer.GetR();
 413
 414     if (i==3 || i==1) {
 415        Double_t rs=0.5*(fgLayers[i+1].GetR() + r);
 416        Double_t d=0.0034, x0=38.6;
 417        if (i==1) {rs=9.; d=0.0097; x0=42;}
 418        if (!fTrackToFollow.PropagateTo(rs,d,x0)) {
 419          //Warning("FollowProlongation","propagation failed !\n");
 420          return;
 421        }
 422     }
 423
 424     //find intersection
 425     Double_t x,y,z;
 426     if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
 427       //Warning("FollowProlongation","failed to estimate track !\n");
 428       return;
 429     }
 430     Double_t phi=TMath::ATan2(y,x);
 431
 432     Int_t idet=layer.FindDetectorIndex(phi,z);
 433     if (idet<0) {
 434       //Warning("FollowProlongation","failed to find a detector !\n");
 435       return;
 436     }
 437
 438     //propagate to the intersection
 439     const AliITSdetector &det=layer.GetDetector(idet);
 440     phi=det.GetPhi();
 441     if (!fTrackToFollow.Propagate(phi,det.GetR())) {
 442       //Warning("FollowProlongation","propagation failed !\n");
 443       return;
 444     }
 445     fTrackToFollow.SetDetectorIndex(idet);
 446
 447     //Select possible prolongations and store the current track estimation
 448     track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
 449     Double_t dz=7*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
 450     Double_t dy=7*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
 451     Double_t road=layer.GetRoad();
 452     if (dz*dy>road*road) {
 453        Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd;
 454        dz=road*scz; dy=road*scy;
 455     }
 456
 457     //Double_t dz=4*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[i]);
 458     if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
 459     if (dz > kMaxRoad) {
 460       //Warning("FollowProlongation","too broad road in Z !\n");
 461       return;
 462     }
 463
 464     if (TMath::Abs(fTrackToFollow.GetZ()-GetZ()) > r+dz) return;
 465
 466     //Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[i]);
 467     if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
 468     if (dy > kMaxRoad) {
 469       //Warning("FollowProlongation","too broad road in Y !\n");
 470       return;
 471     }
 472
 473     fI--;
 474
 475     Double_t zmin=track.GetZ() - dz;
 476     Double_t zmax=track.GetZ() + dz;
 477     Double_t ymin=track.GetY() + r*phi - dy;
 478     Double_t ymax=track.GetY() + r*phi + dy;
 479     if (layer.SelectClusters(zmin,zmax,ymin,ymax)==0)
 480        if (fLayersNotToSkip[fI]) return;
 481
 482     if (!TakeNextProlongation())
 483        if (fLayersNotToSkip[fI]) return;
 484
 485   }
 486
 487   //deal with the best track
 488   Int_t ncl=fTrackToFollow.GetNumberOfClusters();
 489   Int_t nclb=fBestTrack.GetNumberOfClusters();
 490   if (ncl)
 491   if (ncl >= nclb) {
 492      Double_t chi2=fTrackToFollow.GetChi2();
 493      if (chi2/ncl < kChi2PerCluster) {
 494         if (ncl > nclb || chi2 < fBestTrack.GetChi2()) {
 495            ResetBestTrack();
 496         }
 497      }
 498   }
 499
 500 }
 501
 502 Int_t AliITStrackerV2::TakeNextProlongation() {
 503   //--------------------------------------------------------------------
 504   // This function takes another track prolongation
 505   //
 506   //  dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
 507   //--------------------------------------------------------------------
 508   AliITSlayer &layer=fgLayers[fI];
 509   ResetTrackToFollow(fTracks[fI]);
 510
 511   Double_t dz=7*TMath::Sqrt(fTrackToFollow.GetSigmaZ2() + kSigmaZ2[fI]);
 512   Double_t dy=7*TMath::Sqrt(fTrackToFollow.GetSigmaY2() + kSigmaY2[fI]);
 513   Double_t road=layer.GetRoad();
 514   if (dz*dy>road*road) {
 515      Double_t dd=TMath::Sqrt(dz*dy), scz=dz/dd, scy=dy/dd;
 516      dz=road*scz; dy=road*scy;
 517   }
 518
 519   const AliITSRecPoint *c=0; Int_t ci=-1;
 520   const AliITSRecPoint *cc=0; Int_t cci=-1;
 521   Double_t chi2=kMaxChi2;
 522   while ((c=layer.GetNextCluster(ci))!=0) {
 523     Int_t idet=c->GetDetectorIndex();
 524
 525     if (fTrackToFollow.GetDetectorIndex()!=idet) {
 526        const AliITSdetector &det=layer.GetDetector(idet);
 527        ResetTrackToFollow(fTracks[fI]);
 528        if (!fTrackToFollow.Propagate(det.GetPhi(),det.GetR())) {
 529          //Warning("TakeNextProlongation","propagation failed !\n");
 530          continue;
 531        }
 532        fTrackToFollow.SetDetectorIndex(idet);
 533        if (TMath::Abs(fTrackToFollow.GetZ()-GetZ())>layer.GetR()+dz) continue;
 534     }
 535
 536     if (TMath::Abs(fTrackToFollow.GetZ() - c->GetZ()) > dz) continue;
 537     if (TMath::Abs(fTrackToFollow.GetY() - c->GetY()) > dy) continue;
 538
 539     Double_t ch2=fTrackToFollow.GetPredictedChi2(c);
 540     if (ch2 > chi2) continue;
 541     chi2=ch2;
 542     cc=c; cci=ci;
 543     break;
 544   }
 545
 546   if (!cc) return 0;
 547
 548   if (!fTrackToFollow.Update(cc,chi2,(fI<<28)+cci)) {
 549      //Warning("TakeNextProlongation","filtering failed !\n");
 550      return 0;
 551   }
 552
 553   if (fTrackToFollow.GetNumberOfClusters()>1)
 554   if (TMath::Abs(fTrackToFollow.GetD())>4) return 0;
 555
 556   fTrackToFollow.
 557     SetSampledEdx(cc->GetQ(),fTrackToFollow.GetNumberOfClusters()-1); //b.b.
 558
 559   {
 560   Double_t x0;
 561  Double_t d=layer.GetThickness(fTrackToFollow.GetY(),fTrackToFollow.GetZ(),x0);
 562   fTrackToFollow.CorrectForMaterial(d,x0);
 563   }
 564
 565   if (fConstraint[fPass]) {
 566     Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 567     Double_t xyz[]={GetX(),GetY(),GetZ()};
 568     Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 569     fTrackToFollow.Improve(d,xyz,ers);
 570   }
 571
 572   return 1;
 573 }
 574
 575
 576 AliITStrackerV2::AliITSlayer::AliITSlayer() {
 577   //--------------------------------------------------------------------
 578   //default AliITSlayer constructor
 579   //--------------------------------------------------------------------
 580   fR=0.; fPhiOffset=0.; fZOffset=0.;
 581   fNladders=0; fNdetectors=0;
 582   fDetectors=0;
 583
 584   for (Int_t i=0; i<kNsector; i++) fN[i]=0;
 585   fNsel=0;
 586
 587   fRoad=2*fR*TMath::Sqrt(3.14/1.);//assuming that there's only one cluster
 588 }
 589
 590 AliITStrackerV2::AliITSlayer::
 591 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
 592   //--------------------------------------------------------------------
 593   //main AliITSlayer constructor
 594   //--------------------------------------------------------------------
 595   fR=r; fPhiOffset=p; fZOffset=z;
 596   fNladders=nl; fNdetectors=nd;
 597   fDetectors=new AliITSdetector[fNladders*fNdetectors];
 598
 599   for (Int_t i=0; i<kNsector; i++) fN[i]=0;
 600   fNsel=0;
 601
 602   for (Int_t i=0; i<kMaxClusterPerLayer; i++) fClusters[i]=0;
 603
 604   fRoad=2*fR*TMath::Sqrt(3.14/1.);//assuming that there's only one cluster
 605 }
 606
 607 AliITStrackerV2::AliITSlayer::~AliITSlayer() {
 608   //--------------------------------------------------------------------
 609   // AliITSlayer destructor
 610   //--------------------------------------------------------------------
 611   delete[] fDetectors;
 612   ResetClusters();
 613 }
 614
 615 void AliITStrackerV2::AliITSlayer::ResetClusters() {
 616   //--------------------------------------------------------------------
 617   // This function removes loaded clusters
 618   //--------------------------------------------------------------------
 619    for (Int_t s=0; s<kNsector; s++) {
 620        Int_t &n=fN[s];
 621        while (n) {
 622           n--;
 623           delete fClusters[s*kMaxClusterPerSector+n];
 624        }
 625    }
 626 }
 627
 628 void AliITStrackerV2::AliITSlayer::ResetRoad() {
 629   //--------------------------------------------------------------------
 630   // This function calculates the road defined by the cluster density
 631   //--------------------------------------------------------------------
 632   Int_t n=0;
 633   for (Int_t s=0; s<kNsector; s++) {
 634     Int_t i=fN[s];
 635     while (i--)
 636        if (TMath::Abs(fClusters[s*kMaxClusterPerSector+i]->GetZ())<fR) n++;
 637   }
 638   if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
 639 }
 640
 641 Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSRecPoint *c) {
 642   //--------------------------------------------------------------------
 643   // This function inserts a cluster to this layer in increasing
 644   // order of the cluster's fZ
 645   //--------------------------------------------------------------------
 646   Float_t circ=TMath::TwoPi()*fR;
 647   Int_t sec=Int_t(kNsector*c->GetPhiR()/circ);
 648   if (sec>=kNsector) {
 649      ::Error("InsertCluster","Wrong sector !\n");
 650      return 1;
 651   }
 652   Int_t &n=fN[sec];
 653   if (n>=kMaxClusterPerSector) {
 654      ::Error("InsertCluster","Too many clusters !\n");
 655      return 1;
 656   }
 657   if (n==0) fClusters[sec*kMaxClusterPerSector]=c;
 658   else {
 659      Int_t i=FindClusterIndex(c->GetZ(),sec);
 660      Int_t k=n-i+sec*kMaxClusterPerSector;
 661      memmove(fClusters+i+1 ,fClusters+i,k*sizeof(AliITSRecPoint*));
 662      fClusters[i]=c;
 663   }
 664   n++;
 665   return 0;
 666 }
 667
 668 Int_t
 669 AliITStrackerV2::AliITSlayer::FindClusterIndex(Float_t z,Int_t s) const {
 670   //--------------------------------------------------------------------
 671   // For the sector "s", this function returns the index of the first
 672   // with its fZ >= "z".
 673   //--------------------------------------------------------------------
 674   Int_t nc=fN[s];
 675   if (nc==0) return kMaxClusterPerSector*s;
 676
 677   Int_t b=kMaxClusterPerSector*s;
 678   if (z <= fClusters[b]->GetZ()) return b;
 679
 680   Int_t e=b+nc-1;
 681   if (z > fClusters[e]->GetZ()) return e+1;
 682
 683   Int_t m=(b+e)/2;
 684   for (; b<e; m=(b+e)/2) {
 685     if (z > fClusters[m]->GetZ()) b=m+1;
 686     else e=m;
 687   }
 688   return m;
 689 }
 690
 691 Int_t AliITStrackerV2::AliITSlayer::
 692 SelectClusters(Float_t zmin,Float_t zmax,Float_t ymin, Float_t ymax) {
 693   //--------------------------------------------------------------------
 694   // This function selects clusters within the "window"
 695   //--------------------------------------------------------------------
 696     Float_t circ=fR*TMath::TwoPi();
 697
 698     if (ymin>circ) ymin-=circ; else if (ymin<0) ymin+=circ;
 699     if (ymax>circ) ymax-=circ; else if (ymax<0) ymax+=circ;
 700
 701     Int_t i1=Int_t(kNsector*ymin/circ); if (i1==kNsector) i1--;
 702     if (fN[i1]!=0) {
 703        Float_t ym = (ymax<ymin) ? ymax+circ : ymax;
 704        Int_t i=FindClusterIndex(zmin,i1), imax=i1*kMaxClusterPerSector+fN[i1];
 705        for (; i<imax; i++) {
 706            AliITSRecPoint *c=fClusters[i];
 707            if (c->IsUsed()) continue;
 708            if (c->GetZ()>zmax) break;
 709            if (c->GetPhiR()<=ymin) continue;
 710            if (c->GetPhiR()>ym) continue;
 711            fIndex[fNsel++]=i;
 712        }
 713     }
 714
 715     Int_t i2=Int_t(kNsector*ymax/circ); if (i2==kNsector) i2--;
 716     if (i2==i1) return fNsel;
 717
 718     if (fN[i2]!=0) {
 719        Float_t ym = (ymin>ymax) ? ymin-circ : ymin;
 720        Int_t i=FindClusterIndex(zmin,i2), imax=i2*kMaxClusterPerSector+fN[i2];
 721        for (; i<imax; i++) {
 722            AliITSRecPoint *c=fClusters[i];
 723            if (c->IsUsed()) continue;
 724            if (c->GetZ()>zmax) break;
 725            if (c->GetPhiR()<=ym) continue;
 726            if (c->GetPhiR()>ymax) continue;
 727            fIndex[fNsel++]=i;
 728        }
 729     }
 730
 731     return fNsel;
 732 }
 733
 734 const AliITSRecPoint *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
 735   //--------------------------------------------------------------------
 736   // This function returns clusters within the "window"
 737   //--------------------------------------------------------------------
 738   AliITSRecPoint *c=0;
 739   ci=-1;
 740   if (fNsel) {
 741      fNsel--;
 742      ci=fIndex[fNsel];
 743      c=fClusters[ci];
 744   }
 745   return c;
 746 }
 747
 748 Int_t AliITStrackerV2::AliITSlayer::GetNumberOfClusters() const {
 749   Int_t n=0;
 750   for (Int_t s=0; s<kNsector; s++) n+=fN[s];
 751   return n;
 752 }
 753
 754 Int_t
 755 AliITStrackerV2::AliITSlayer::FindDetectorIndex(Double_t phi,Double_t z)const {
 756   //--------------------------------------------------------------------
 757   //This function finds the detector crossed by the track
 758   //--------------------------------------------------------------------
 759   Double_t dphi=-(phi-fPhiOffset);
 760   if      (dphi <  0) dphi += 2*TMath::Pi();
 761   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
 762   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
 763   if (np>=fNladders) np-=fNladders;
 764   if (np<0)          np+=fNladders;
 765
 766   Double_t dz=fZOffset-z;
 767   Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
 768   if (nz>=fNdetectors) return -1;
 769   if (nz<0)            return -1;
 770
 771   return np*fNdetectors + nz;
 772 }
 773
 774 Double_t
 775 AliITStrackerV2::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
 776 const {
 777   //--------------------------------------------------------------------
 778   //This function returns the layer thickness at this point (units X0)
 779   //--------------------------------------------------------------------
 780   Double_t d=0.0085;
 781   x0=21.82;
 782
 783   if (43<fR&&fR<45) { //SSD2
 784      Double_t dd=0.0034;
 785      d=dd;
 786      if (TMath::Abs(y-0.00)>3.40) d+=dd;
 787      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
 788      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
 789      for (Int_t i=0; i<12; i++) {
 790        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
 791           if (TMath::Abs(y-0.00)>3.40) d+=dd;
 792           d+=0.0034;
 793           break;
 794        }
 795        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
 796           if (TMath::Abs(y-0.00)>3.40) d+=dd;
 797           d+=0.0034;
 798           break;
 799        }
 800        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
 801        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
 802      }
 803   } else
 804   if (37<fR&&fR<41) { //SSD1
 805      Double_t dd=0.0034;
 806      d=dd;
 807      if (TMath::Abs(y-0.00)>3.40) d+=dd;
 808      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
 809      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
 810      for (Int_t i=0; i<11; i++) {
 811        if (TMath::Abs(z-3.9*i)<0.15) {
 812           if (TMath::Abs(y-0.00)>3.40) d+=dd;
 813           d+=dd;
 814           break;
 815        }
 816        if (TMath::Abs(z+3.9*i)<0.15) {
 817           if (TMath::Abs(y-0.00)>3.40) d+=dd;
 818           d+=dd;
 819           break;
 820        }
 821        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
 822        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
 823      }
 824   } else
 825   if (13<fR&&fR<26) { //SDD
 826      Double_t dd=0.0033;
 827      d=dd;
 828      if (TMath::Abs(y-0.00)>3.30) d+=dd;
 829
 830      if (TMath::Abs(y-1.80)<0.55) {
 831         d+=0.016;
 832         for (Int_t j=0; j<20; j++) {
 833           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
 834           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
 835         }
 836      }
 837      if (TMath::Abs(y+1.80)<0.55) {
 838         d+=0.016;
 839         for (Int_t j=0; j<20; j++) {
 840           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
 841           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
 842         }
 843      }
 844
 845      for (Int_t i=0; i<4; i++) {
 846        if (TMath::Abs(z-7.3*i)<0.60) {
 847           d+=dd;
 848           if (TMath::Abs(y-0.00)>3.30) d+=dd;
 849           break;
 850        }
 851        if (TMath::Abs(z+7.3*i)<0.60) {
 852           d+=dd;
 853           if (TMath::Abs(y-0.00)>3.30) d+=dd;
 854           break;
 855        }
 856      }
 857   } else
 858   if (6<fR&&fR<8) {   //SPD2
 859      Double_t dd=0.0063; x0=21.5;
 860      d=dd;
 861      if (TMath::Abs(y-3.08)>0.5) d+=dd;
 862      //if (TMath::Abs(y-3.08)>0.45) d+=dd;
 863      if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
 864   } else
 865   if (3<fR&&fR<5) {   //SPD1
 866      Double_t dd=0.0063; x0=21.5;
 867      d=dd;
 868      if (TMath::Abs(y+0.21)>0.6) d+=dd;
 869      //if (TMath::Abs(y+0.21)>0.45) d+=dd;
 870      if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
 871   }
 872
 873   return d;
 874 }
 875
 876 Double_t AliITStrackerV2::GetEffectiveThickness(Double_t y,Double_t z) const
 877 {
 878   //--------------------------------------------------------------------
 879   //Returns the thickness between the current layer and the vertex (units X0)
 880   //--------------------------------------------------------------------
 881   Double_t d=0.0028*3*3; //beam pipe
 882   Double_t x0=0;
 883
 884   Double_t xn=fgLayers[fI].GetR();
 885   for (Int_t i=0; i<fI; i++) {
 886     Double_t xi=fgLayers[i].GetR();
 887     d+=fgLayers[i].GetThickness(y,z,x0)*xi*xi;
 888   }
 889
 890   if (fI>1) {
 891     Double_t xi=9.;
 892     d+=0.0097*xi*xi;
 893   }
 894
 895   if (fI>3) {
 896     Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR());
 897     d+=0.0034*xi*xi;
 898   }
 899
 900   return d/(xn*xn);
 901 }
 902
 903 Bool_t
 904 AliITStrackerV2::RefitAt(Double_t xx,AliITStrackV2 *t,const AliITStrackV2 *c) {
 905   //--------------------------------------------------------------------
 906   // This function refits the track "t" at the position "x" using
 907   // the clusters from "c"
 908   //--------------------------------------------------------------------
 909   Int_t index[kMaxLayer];
 910   Int_t k;
 911   for (k=0; k<kMaxLayer; k++) index[k]=-1;
 912   Int_t nc=c->GetNumberOfClusters();
 913   for (k=0; k<nc; k++) {
 914     Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
 915     index[nl]=idx;
 916   }
 917
 918   Int_t from, to, step;
 919   if (xx > t->GetX()) {
 920       from=0; to=kMaxLayer;
 921       step=+1;
 922   } else {
 923       from=kMaxLayer-1; to=-1;
 924       step=-1;
 925   }
 926
 927   for (Int_t i=from; i != to; i += step) {
 928      AliITSlayer &layer=fgLayers[i];
 929      Double_t r=layer.GetR();
 930
 931      {
 932      Double_t hI=i-0.5*step;
 933      if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
 934         Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
 935         Double_t d=0.0034, x0=38.6;
 936         if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
 937         if (!t->PropagateTo(rs,-step*d,x0)) {
 938           return kFALSE;
 939         }
 940      }
 941      }
 942
 943      // remember old position [SR, GSI 18.02.2003]
 944      Double_t oldX=0., oldY=0., oldZ=0.;
 945      if (t->IsStartedTimeIntegral() && step==1) {
 946         t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
 947      }
 948      //
 949
 950      Double_t x,y,z;
 951      if (!t->GetGlobalXYZat(r,x,y,z)) {
 952        return kFALSE;
 953      }
 954      Double_t phi=TMath::ATan2(y,x);
 955      Int_t idet=layer.FindDetectorIndex(phi,z);
 956      if (idet<0) {
 957        return kFALSE;
 958      }
 959      const AliITSdetector &det=layer.GetDetector(idet);
 960      phi=det.GetPhi();
 961      if (!t->Propagate(phi,det.GetR())) {
 962        return kFALSE;
 963      }
 964      t->SetDetectorIndex(idet);
 965
 966      const AliITSRecPoint *cl=0;
 967      Double_t maxchi2=kMaxChi2;
 968
 969      Int_t idx=index[i];
 970      if (idx>0) {
 971         const AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(idx);
 972         if (idet != c->GetDetectorIndex()) {
 973            idet=c->GetDetectorIndex();
 974            const AliITSdetector &det=layer.GetDetector(idet);
 975            if (!t->Propagate(det.GetPhi(),det.GetR())) {
 976              return kFALSE;
 977            }
 978            t->SetDetectorIndex(idet);
 979         }
 980         Double_t chi2=t->GetPredictedChi2(c);
 981         if (chi2<maxchi2) {
 982           cl=c;
 983           maxchi2=chi2;
 984         } else {
 985           return kFALSE;
 986         }
 987      }
 988      /*
 989      if (cl==0)
 990      if (t->GetNumberOfClusters()>2) {
 991         Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
 992         Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
 993         Double_t zmin=t->GetZ() - dz;
 994         Double_t zmax=t->GetZ() + dz;
 995         Double_t ymin=t->GetY() + phi*r - dy;
 996         Double_t ymax=t->GetY() + phi*r + dy;
 997         layer.SelectClusters(zmin,zmax,ymin,ymax);
 998
 999         const AliITSRecPoint *c=0; Int_t ci=-1;
1000         while ((c=layer.GetNextCluster(ci))!=0) {
1001            if (idet != c->GetDetectorIndex()) continue;
1002            Double_t chi2=t->GetPredictedChi2(c);
1003            if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1004         }
1005      }
1006      */
1007      if (cl) {
1008        if (!t->Update(cl,maxchi2,idx)) {
1009           return kFALSE;
1010        }
1011        t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1012      }
1013
1014      {
1015      Double_t x0;
1016      Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1017      t->CorrectForMaterial(-step*d,x0);
1018      }
1019
1020      // track time update [SR, GSI 17.02.2003]
1021      if (t->IsStartedTimeIntegral() && step==1) {
1022         Double_t newX, newY, newZ;
1023         t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1024         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1025                        (oldZ-newZ)*(oldZ-newZ);
1026         t->AddTimeStep(TMath::Sqrt(dL2));
1027      }
1028      //
1029
1030   }
1031
1032   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1033   return kTRUE;
1034 }
1035
1036 void AliITStrackerV2::UseClusters(const AliKalmanTrack *t, Int_t from) const {
1037   //--------------------------------------------------------------------
1038   // This function marks clusters assigned to the track
1039   //--------------------------------------------------------------------
1040   AliTracker::UseClusters(t,from);
1041
1042   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
1043   if (c && c->GetSigmaZ2()>0.1) c->UnUse();
1044   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
1045   if (c && c->GetSigmaZ2()>0.1) c->UnUse();
1046
1047 }