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