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