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