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 //-------------------------------------------------------------------------
  21 #include <TFile.h>
  22 #include <TTree.h>
  23 #include <TRandom.h>
  24 #include <iostream.h>
  25
  26 #include "AliITSgeom.h"
  27 #include "AliITSRecPoint.h"
  28 #include "AliTPCtrack.h"
  29 #include "AliITSclusterV2.h"
  30 #include "AliITStrackerV2.h"
  31
  32 //#define DEBUG
  33
  34 #ifdef DEBUG
  35 Int_t LAB=236;
  36 #endif
  37
  38 extern TRandom *gRandom;
  39
  40 AliITStrackerV2::AliITSlayer AliITStrackerV2::fLayers[kMaxLayer]; // ITS layers
  41
  42 AliITStrackerV2::
  43 AliITStrackerV2(const AliITSgeom *geom) throw (const Char_t *) {
  44   //--------------------------------------------------------------------
  45   //This is the AliITStracker constructor
  46   //It also reads clusters from a root file
  47   //--------------------------------------------------------------------
  48   fYV=fZV=0.;
  49
  50   AliITSgeom *g=(AliITSgeom*)geom;
  51
  52   Float_t x,y,z;
  53   Int_t i;
  54   for (i=1; i<kMaxLayer+1; i++) {
  55     Int_t nlad=g->GetNladders(i);
  56     Int_t ndet=g->GetNdetectors(i);
  57
  58     g->GetTrans(i,1,1,x,y,z);
  59     Double_t r=TMath::Sqrt(x*x + y*y);
  60     Double_t poff=TMath::ATan2(y,x);
  61     Double_t zoff=z;
  62
  63     g->GetTrans(i,1,2,x,y,z);
  64     r += TMath::Sqrt(x*x + y*y);
  65     g->GetTrans(i,2,1,x,y,z);
  66     r += TMath::Sqrt(x*x + y*y);
  67     g->GetTrans(i,2,2,x,y,z);
  68     r += TMath::Sqrt(x*x + y*y);
  69     r*=0.25;
  70
  71     new (fLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
  72
  73     for (Int_t j=1; j<nlad+1; j++) {
  74       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
  75         Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift);
  76         Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
  77
  78         Double_t r =-x*rot[1] + y*rot[0];         if (i==1) r=-r;
  79         Double_t phi=TMath::ATan2(rot[1],rot[0]); if (i==1) phi-=3.1415927;
  80         phi+=0.5*TMath::Pi();
  81         AliITSdetector &det=fLayers[i-1].GetDetector((j-1)*ndet + k-1);
  82
  83 if (phi<0) phi += 2*TMath::Pi();
  84
  85         new(&det) AliITSdetector(r,phi);
  86       }
  87     }
  88
  89   }
  90
  91   try {
  92      //Read clusters
  93      TTree *cTree=(TTree*)gDirectory->Get("cTree");
  94      if (!cTree) throw
  95         ("AliITStrackerV2::AliITStrackerV2 can't get cTree !\n");
  96
  97      TBranch *branch=cTree->GetBranch("Clusters");
  98      if (!branch) throw
  99         ("AliITStrackerV2::AliITStrackerV2 can't get Clusters branch !\n");
 100
 101      TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
 102      branch->SetAddress(&clusters);
 103
 104      Int_t nentr=(Int_t)cTree->GetEntries();
 105      for (i=0; i<nentr; i++) {
 106        if (!cTree->GetEvent(i)) continue;
 107        //Int_t lay,lad,det; g->GetModuleId(i-1,lay,lad,det);
 108        Int_t lay,lad,det; g->GetModuleId(i,lay,lad,det);
 109        Int_t ncl=clusters->GetEntriesFast();
 110        while (ncl--) {
 111          AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
 112
 113 #ifdef DEBUG
 114 if (c->GetLabel(2)!=LAB)
 115 if (c->GetLabel(1)!=LAB)
 116 if (c->GetLabel(0)!=LAB) continue;
 117 cout<<lay-1<<' '<<lad-1<<' '<<det-1<<' '<<c->GetY()<<' '<<c->GetZ()<<endl;
 118 #endif
 119
 120          AliITSdetector &det=fLayers[lay-1].GetDetector(c->GetDetectorIndex());
 121          Double_t r=det.GetR(); r=TMath::Sqrt(r*r + c->GetY()*c->GetY());
 122          if (r > TMath::Abs(c->GetZ())-0.5)
 123             fLayers[lay-1].InsertCluster(new AliITSclusterV2(*c));
 124        }
 125        clusters->Delete();
 126      }
 127   }
 128   catch (const Char_t *msg) {
 129     cerr<<msg<<endl;
 130     throw;
 131   }
 132
 133   fI=kMaxLayer;
 134 }
 135
 136 #ifdef DEBUG
 137 static Int_t lbl;
 138 #endif
 139
 140 Int_t AliITStrackerV2::Clusters2Tracks(const TFile *inp, TFile *out) {
 141   //--------------------------------------------------------------------
 142   //This functions reconstructs ITS tracks
 143   //--------------------------------------------------------------------
 144   TFile *in=(TFile*)inp;
 145   TDirectory *savedir=gDirectory;
 146
 147   if (!in->IsOpen()) {
 148      cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
 149      cerr<<"file with TPC tracks is not open !\n";
 150      return 1;
 151   }
 152
 153   if (!out->IsOpen()) {
 154      cerr<<"AliITStrackerV2::Clusters2Tracks(): ";
 155      cerr<<"file for ITS tracks is not open !\n";
 156      return 2;
 157   }
 158
 159   in->cd();
 160   TTree *tpcTree=(TTree*)in->Get("TPCf");
 161   if (!tpcTree) {
 162      cerr<<"AliITStrackerV2::Clusters2Tracks() ";
 163      cerr<<"can't get a tree with TPC tracks !\n";
 164      return 3;
 165   }
 166
 167   AliTPCtrack *itrack=new AliTPCtrack;
 168   tpcTree->SetBranchAddress("tracks",&itrack);
 169
 170   out->cd();
 171   TTree itsTree("ITSf","Tree with ITS tracks");
 172   AliITStrackV2 *otrack=&fBestTrack;
 173   itsTree.Branch("tracks","AliITStrackV2",&otrack,32000,0);
 174
 175   Int_t ntrk=0;
 176
 177   Int_t nentr=(Int_t)tpcTree->GetEntries();
 178   for (Int_t i=0; i<nentr; i++) {
 179
 180     if (!tpcTree->GetEvent(i)) continue;
 181     Int_t tpcLabel=itrack->GetLabel(); //save the TPC track label
 182
 183 #ifdef DEBUG
 184 lbl=tpcLabel;
 185 if (TMath::Abs(tpcLabel)!=LAB) continue;
 186 cout<<tpcLabel<<" *****************\n";
 187 #endif
 188
 189     try {
 190       ResetTrackToFollow(AliITStrackV2(*itrack));
 191     } catch (const Char_t *msg) {
 192       cerr<<msg<<endl;
 193       continue;
 194     }
 195     ResetBestTrack();
 196
 197     fYV=gRandom->Gaus(0.,kSigmaYV); fZV=gRandom->Gaus(0.,kSigmaZV);
 198
 199     Double_t r2=fTrackToFollow.GetX()*fTrackToFollow.GetX();
 200     Double_t x0=0.082/21.82*2.33*(45*45+40*40)/r2+2.000/41.28*0.03*75*75/r2;
 201     fTrackToFollow.Improve(x0,fYV,fZV);
 202
 203     Double_t xk=80.;
 204     fTrackToFollow.PropagateTo(xk,0.,0.); //Ne if it's still there
 205
 206     xk-=0.005;
 207     fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
 208     xk-=0.02;
 209     fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);   //Kevlar
 210        xk-=2.0;
 211        fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
 212     xk-=0.02;
 213     fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);   //Kevlar
 214     xk-=0.005;
 215     fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
 216
 217     xk=61.;
 218     fTrackToFollow.PropagateTo(xk,0.,0.); //C02
 219
 220     xk -=0.005;
 221     fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7);    //Al
 222     xk -=0.005;
 223     fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71);  //Tedlar
 224     xk -=0.02;
 225     fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);    //Kevlar
 226        xk -=0.5;
 227        fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029);  //Nomex
 228     xk -=0.02;
 229     fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);    //Kevlar
 230     xk -=0.005;
 231     fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71);  //Tedlar
 232     xk -=0.005;
 233     fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7);    //Al
 234
 235
 236     for (FollowProlongation(); fI<kMaxLayer; fI++) {
 237        while (TakeNextProlongation()) FollowProlongation();
 238     }
 239
 240 #ifdef DEBUG
 241 cout<<fBestTrack.GetNumberOfClusters()<<" number of clusters\n\n";
 242 #endif
 243
 244     if (fBestTrack.GetNumberOfClusters() >= kMaxLayer-kLayersToSkip) {
 245        cerr << ++ntrk << "                \r";
 246        fBestTrack.SetLabel(tpcLabel);
 247        UseClusters(&fBestTrack);
 248        itsTree.Fill();
 249     }
 250
 251   }
 252
 253   itsTree.Write();
 254   savedir->cd();
 255   cerr<<"Number of TPC tracks: "<<nentr<<endl;
 256   cerr<<"Number of prolonged tracks: "<<ntrk<<endl;
 257
 258   delete itrack;
 259
 260   return 0;
 261 }
 262
 263 Int_t AliITStrackerV2::PropagateBack(const TFile *inp, TFile *out) {
 264   //--------------------------------------------------------------------
 265   //This functions propagates reconstructed ITS tracks back
 266   //--------------------------------------------------------------------
 267   TFile *in=(TFile*)inp;
 268   TDirectory *savedir=gDirectory;
 269
 270   if (!in->IsOpen()) {
 271      cerr<<"AliITStrackerV2::PropagateBack(): ";
 272      cerr<<"file with ITS tracks is not open !\n";
 273      return 1;
 274   }
 275
 276   if (!out->IsOpen()) {
 277      cerr<<"AliITStrackerV2::PropagateBack(): ";
 278      cerr<<"file for back propagated ITS tracks is not open !\n";
 279      return 2;
 280   }
 281
 282   in->cd();
 283   TTree *itsTree=(TTree*)in->Get("ITSf");
 284   if (!itsTree) {
 285      cerr<<"AliITStrackerV2::PropagateBack() ";
 286      cerr<<"can't get a tree with ITS tracks !\n";
 287      return 3;
 288   }
 289   AliITStrackV2 *itrack=new AliITStrackV2;
 290   itsTree->SetBranchAddress("tracks",&itrack);
 291
 292   out->cd();
 293   TTree backTree("ITSb","Tree with back propagated ITS tracks");
 294   AliTPCtrack *otrack=0;
 295   backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
 296
 297   Int_t ntrk=0;
 298
 299   Int_t nentr=(Int_t)itsTree->GetEntries();
 300   for (Int_t i=0; i<nentr; i++) {
 301     itsTree->GetEvent(i);
 302     ResetTrackToFollow(*itrack);
 303     fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
 304     Int_t itsLabel=fTrackToFollow.GetLabel(); //save the ITS track label
 305
 306 #ifdef DEBUG
 307 if (TMath::Abs(itsLabel)!=LAB) continue;
 308 cout<<itsLabel<<" *****************\n";
 309 #endif
 310
 311     try {
 312        Int_t nc=itrack->GetNumberOfClusters();
 313 #ifdef DEBUG
 314 for (Int_t k=0; k<nc; k++) {
 315     Int_t index=itrack->GetClusterIndex(k);
 316     AliITSclusterV2 *c=(AliITSclusterV2*)GetCluster(index);
 317     cout<<c->GetLabel(0)<<' '<<c->GetLabel(1)<<' '<<c->GetLabel(2)<<endl;
 318 }
 319 #endif
 320        Int_t idx=0, l=0;
 321        const  AliITSclusterV2 *c=0;
 322        if (nc--) {
 323           idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
 324           c=(AliITSclusterV2*)GetCluster(idx);
 325        }
 326        for (fI=0; fI<kMaxLayer; fI++) {
 327          AliITSlayer &layer=fLayers[fI];
 328          Double_t r=layer.GetR();
 329          if (fI==2 || fI==4) {
 330             Double_t rs=0.5*(fLayers[fI-1].GetR() + r);
 331             Double_t ds=0.034; if (fI==4) ds=0.039;
 332             Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
 333             if (!fTrackToFollow.PropagateTo(rs,1*dx0r,dr)) throw "";
 334          }
 335
 336          Double_t x,y,z;
 337          if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z))
 338             throw "AliITStrackerV2::PropagateBack: failed to estimate track !";
 339          Double_t phi=TMath::ATan2(y,x);
 340          Double_t d=layer.GetThickness(phi,z);
 341          Int_t idet=layer.FindDetectorIndex(phi,z);
 342          if (idet<0)
 343          throw "AliITStrackerV2::PropagateBack: failed to find a detector !\n";
 344          const AliITSdetector &det=layer.GetDetector(idet);
 345          r=det.GetR(); phi=det.GetPhi();
 346          if (!fTrackToFollow.Propagate(phi,r,d/21.82*2.33,d*2.33)) throw "";
 347
 348          const AliITSclusterV2 *cl=0;
 349          Int_t index=0;
 350          Double_t maxchi2=kMaxChi2;
 351
 352          if (l==fI) {
 353            if (idet != c->GetDetectorIndex()) {
 354              idet=c->GetDetectorIndex();
 355              const AliITSdetector &det=layer.GetDetector(idet);
 356              r=det.GetR(); phi=det.GetPhi();
 357              if (!fTrackToFollow.Propagate(phi,r,0.,0.)) throw "";
 358            }
 359            Double_t chi2=fTrackToFollow.GetPredictedChi2(c);
 360            if (chi2<kMaxChi2) {
 361               cl=c; maxchi2=chi2; index=idx;
 362            }
 363            if (nc--) {
 364               idx=itrack->GetClusterIndex(nc); l=(idx&0xf0000000)>>28;
 365               c=(AliITSclusterV2*)GetCluster(idx);
 366            }
 367          }
 368
 369          if (fTrackToFollow.GetNumberOfClusters()>2) {
 370            Double_t dz=3*TMath::Sqrt(fTrackToFollow.GetSigmaZ2()+kSigmaZ2[fI]);
 371            Double_t dy=3*TMath::Sqrt(fTrackToFollow.GetSigmaY2()+kSigmaY2[fI]);
 372            Double_t zmin=fTrackToFollow.GetZ() - dz;
 373            Double_t zmax=fTrackToFollow.GetZ() + dz;
 374            Double_t ymin=fTrackToFollow.GetY() + phi*r - dy;
 375            Double_t ymax=fTrackToFollow.GetY() + phi*r + dy;
 376            layer.SelectClusters(zmin,zmax,ymin,ymax);
 377
 378            const AliITSclusterV2 *cc=0; Int_t ci;
 379            while ((cc=layer.GetNextCluster(ci))!=0) {
 380               if (idet != cc->GetDetectorIndex()) continue;
 381               Double_t chi2=fTrackToFollow.GetPredictedChi2(cc);
 382               if (chi2<maxchi2) {
 383                  cl=cc; index=(fI<<28)+ci; maxchi2=chi2;
 384               }
 385            }
 386          }
 387
 388          if (cl) {
 389             if (!fTrackToFollow.Update(cl,maxchi2,index))
 390               cerr<<"AliITStrackerV2::PropagateBack: filtering failed !\n";
 391               continue;
 392          }
 393        }
 394
 395        Double_t xk=61.;
 396        fTrackToFollow.PropagateTo(xk,0.,0.); //Air
 397
 398        xk +=0.005;
 399        fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7);    //Al
 400        xk +=0.005;
 401        fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71);  //Tedlar
 402        xk +=0.02;
 403        fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);    //Kevlar
 404           xk +=0.5;
 405           fTrackToFollow.PropagateTo(xk, 0.5/41.28*.029, 0.5*0.029);  //Nomex
 406        xk +=0.02;
 407        fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);    //Kevlar
 408        xk +=0.005;
 409        fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71);  //Tedlar
 410        xk +=0.005;
 411        fTrackToFollow.PropagateTo(xk, 0.005/24.01*2.7, 0.005*2.7);    //Al
 412
 413        xk=80.;
 414        fTrackToFollow.PropagateTo(xk,0.,0.); //CO2
 415
 416        xk+=0.005;
 417        fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
 418        xk+=0.02;
 419        fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);   //Kevlar
 420           xk+=2.0;
 421           fTrackToFollow.PropagateTo(xk, 2.0/41.28*0.029, 2.0*0.029);//Nomex
 422        xk+=0.02;
 423        fTrackToFollow.PropagateTo(xk, 0.02/44.86*1.45, 0.02*1.45);   //Kevlar
 424        xk+=0.005;
 425        fTrackToFollow.PropagateTo(xk, 0.005/44.77*1.71, 0.005*1.71); //Tedlar
 426
 427        fTrackToFollow.SetLabel(itsLabel);
 428        otrack=new AliTPCtrack(fTrackToFollow,fTrackToFollow.GetAlpha());
 429        backTree.Fill(); delete otrack;
 430        UseClusters(&fTrackToFollow);
 431        cerr << ++ntrk << "                \r";
 432     }
 433     catch (const Char_t *msg) {
 434        cerr<<msg<<endl;
 435     }
 436   }
 437
 438   backTree.Write();
 439   savedir->cd();
 440   cerr<<"Number of ITS tracks: "<<nentr<<endl;
 441   cerr<<"Number of back propagated ITS tracks: "<<ntrk<<endl;
 442
 443   delete itrack;
 444
 445   return 0;
 446 }
 447
 448 AliCluster *AliITStrackerV2::GetCluster(Int_t index) const {
 449   //--------------------------------------------------------------------
 450   //       Return pointer to a given cluster
 451   //--------------------------------------------------------------------
 452   Int_t l=(index & 0xf0000000) >> 28;
 453   Int_t c=(index & 0x0fffffff) >> 00;
 454   return fLayers[l].GetCluster(c);
 455 }
 456
 457
 458 void AliITStrackerV2::FollowProlongation() {
 459   //--------------------------------------------------------------------
 460   //This function finds a track prolongation
 461   //--------------------------------------------------------------------
 462   Int_t tryAgain=kLayersToSkip;
 463
 464   while (fI) {
 465     fI--;
 466 #ifdef DEBUG
 467 cout<<fI<<' ';
 468 #endif
 469     AliITSlayer &layer=fLayers[fI];
 470     AliITStrackV2 &track=fTracks[fI];
 471
 472     Double_t r=layer.GetR();
 473     if (fI==3 || fI==1) {
 474       Double_t rs=0.5*(fLayers[fI+1].GetR() + r);
 475       Double_t ds=0.034; if (fI==3) ds=0.039;
 476       Double_t dx0r=ds/21.82*2.33, dr=ds*2.33;
 477       fTrackToFollow.Propagate(fTrackToFollow.GetAlpha(),rs,1*dx0r,dr);
 478     }
 479
 480     //find intersection
 481     Double_t x,y,z;
 482     if (!fTrackToFollow.GetGlobalXYZat(r,x,y,z)) {
 483      cerr<<"AliITStrackerV2::FollowProlongation: failed to estimate track !\n";
 484      break;
 485     }
 486     Double_t phi=TMath::ATan2(y,x);
 487     Double_t d=layer.GetThickness(phi,z);
 488     Int_t idet=layer.FindDetectorIndex(phi,z);
 489     if (idet<0) {
 490     cerr<<"AliITStrackerV2::FollowProlongation: failed to find a detector !\n";
 491       break;
 492     }
 493
 494     //propagate to the intersection
 495     const AliITSdetector &det=layer.GetDetector(idet);
 496     phi=det.GetPhi();
 497     if (!fTrackToFollow.Propagate(phi,det.GetR(),1*d/21.82*2.33,d*2.33)) {
 498       cerr<<"AliITStrackerV2::FollowProlongation: propagation failed !\n";
 499       break;
 500     }
 501     fTrackToFollow.SetDetectorIndex(idet);
 502
 503     //Select possible prolongations and store the current track estimation
 504     track.~AliITStrackV2(); new(&track) AliITStrackV2(fTrackToFollow);
 505     Double_t dz=3*TMath::Sqrt(track.GetSigmaZ2() + kSigmaZ2[fI]);
 506     if (dz < 0.5*TMath::Abs(track.GetTgl())) dz=0.5*TMath::Abs(track.GetTgl());
 507     if (dz > kMaxRoad/4) {
 508       //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Z !\n";
 509       break;
 510     }
 511     Double_t dy=4*TMath::Sqrt(track.GetSigmaY2() + kSigmaY2[fI]);
 512     if (dy < 0.5*TMath::Abs(track.GetSnp())) dy=0.5*TMath::Abs(track.GetSnp());
 513     if (dy > kMaxRoad/4) {
 514       //cerr<<"AliITStrackerV2::FollowProlongation: too broad road in Y !\n";
 515       break;
 516     }
 517
 518     Double_t zmin=track.GetZ() - dz;
 519     Double_t zmax=track.GetZ() + dz;
 520     Double_t ymin=track.GetY() + r*phi - dy;
 521     Double_t ymax=track.GetY() + r*phi + dy;
 522     layer.SelectClusters(zmin,zmax,ymin,ymax);
 523
 524     //take another prolongation
 525     if (!TakeNextProlongation()) if (!tryAgain--) break;
 526     tryAgain=kLayersToSkip;
 527
 528   }
 529
 530   //deal with the best track
 531   Int_t ncl=fTrackToFollow.GetNumberOfClusters();
 532   Int_t nclb=fBestTrack.GetNumberOfClusters();
 533   if (ncl)
 534   if (ncl >= nclb) {
 535      Double_t chi2=fTrackToFollow.GetChi2();
 536      if (chi2/ncl < kChi2PerCluster) {
 537         if (ncl > nclb || chi2 < fBestTrack.GetChi2()) {
 538            ResetBestTrack();
 539         }
 540      }
 541   }
 542
 543   if (fI) fI++;
 544 }
 545
 546
 547 Int_t AliITStrackerV2::TakeNextProlongation() {
 548   //--------------------------------------------------------------------
 549   //This function takes another track prolongation
 550   //--------------------------------------------------------------------
 551   //Double_t m[20];
 552   Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 553
 554   AliITSlayer &layer=fLayers[fI];
 555   AliITStrackV2 &t=fTracks[fI];
 556
 557   Double_t dz=4*TMath::Sqrt(t.GetSigmaZ2() + kSigmaZ2[fI]);
 558   Double_t dy=4*TMath::Sqrt(t.GetSigmaY2() + kSigmaY2[fI]);
 559
 560   const AliITSclusterV2 *c=0; Int_t ci=-1;
 561   Double_t chi2=12345.;
 562   while ((c=layer.GetNextCluster(ci))!=0) {
 563
 564 #ifdef DEBUG
 565 //if (fI==0)
 566 //if (c->GetLabel(0)!=TMath::Abs(lbl)) continue;
 567 #endif
 568
 569     Int_t idet=c->GetDetectorIndex();
 570
 571     if (t.GetDetectorIndex()!=idet) {
 572        const AliITSdetector &det=layer.GetDetector(idet);
 573        if (!t.Propagate(det.GetPhi(),det.GetR(),0.,0.)) {
 574          cerr<<"AliITStrackerV2::TakeNextProlongation: propagation failed !\n";
 575          continue;
 576        }
 577        t.SetDetectorIndex(idet);
 578
 579 #ifdef DEBUG
 580 cout<<fI<<" change detector !\n";
 581 #endif
 582
 583     }
 584
 585     if (TMath::Abs(t.GetZ() - c->GetZ()) > dz) continue;
 586     if (TMath::Abs(t.GetY() - c->GetY()) > dy) continue;
 587
 588     //m[0]=fYV; m[1]=fZV;
 589     //chi2=t.GetPredictedChi2(c,m,d/21.82*2.33);
 590     chi2=t.GetPredictedChi2(c);
 591
 592     if (chi2<kMaxChi2) break;
 593   }
 594
 595 #ifdef DEBUG
 596 cout<<fI<<" chi2="<<chi2<<' '<<t.GetY()<<' '<<t.GetZ()<<' '<<dy<<' '<<dz<<endl;
 597 #endif
 598
 599   if (chi2>=kMaxChi2) return 0;
 600   if (!c) return 0;
 601
 602   ResetTrackToFollow(t);
 603
 604   //if (!fTrackToFollow.Update(m,chi2,(fI<<28)+ci)) {
 605   if (!fTrackToFollow.Update(c,chi2,(fI<<28)+ci)) {
 606      cerr<<"AliITStrackerV2::TakeNextProlongation: filtering failed !\n";
 607      return 0;
 608   }
 609   fTrackToFollow.Improve(d/21.82*2.33,fYV,fZV);
 610
 611 #ifdef DEBUG
 612 cout<<"accepted lab="<<c->GetLabel(0)<<' '
 613     <<fTrackToFollow.GetNumberOfClusters()<<' '
 614     <<fTrackToFollow.GetY()<<' '<<fTrackToFollow.GetZ()<<endl<<endl;
 615 #endif
 616
 617   return 1;
 618 }
 619
 620
 621
 622
 623 AliITStrackerV2::AliITSlayer::AliITSlayer() {
 624   //--------------------------------------------------------------------
 625   //default AliITSlayer constructor
 626   //--------------------------------------------------------------------
 627   fN=0;
 628   fDetectors=0;
 629 }
 630
 631 AliITStrackerV2::AliITSlayer::
 632 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
 633   //--------------------------------------------------------------------
 634   //main AliITSlayer constructor
 635   //--------------------------------------------------------------------
 636   fR=r; fPhiOffset=p; fZOffset=z;
 637   fNladders=nl; fNdetectors=nd;
 638   fDetectors=new AliITSdetector[fNladders*fNdetectors];
 639
 640   fN=0;
 641   fI=0;
 642 }
 643
 644 AliITStrackerV2::AliITSlayer::~AliITSlayer() {
 645   //--------------------------------------------------------------------
 646   // AliITSlayer destructor
 647   //--------------------------------------------------------------------
 648   delete[] fDetectors;
 649   for (Int_t i=0; i<fN; i++) delete fClusters[i];
 650 }
 651
 652 Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
 653   //--------------------------------------------------------------------
 654   //This function adds a cluster to this layer
 655   //--------------------------------------------------------------------
 656   if (fN==kMaxClusterPerLayer) {
 657  cerr<<"AliITStrackerV2::AliITSlayer::InsertCluster(): Too many clusters !\n";
 658    return 1;
 659   }
 660
 661   if (fN==0) {fClusters[fN++]=c; return 0;}
 662   Int_t i=FindClusterIndex(c->GetZ());
 663   memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliITSclusterV2*));
 664   fClusters[i]=c; fN++;
 665
 666   return 0;
 667 }
 668
 669 Int_t AliITStrackerV2::AliITSlayer::FindClusterIndex(Double_t z) const {
 670   //--------------------------------------------------------------------
 671   // This function returns the index of the nearest cluster
 672   //--------------------------------------------------------------------
 673   if (fN==0) return 0;
 674   if (z <= fClusters[0]->GetZ()) return 0;
 675   if (z > fClusters[fN-1]->GetZ()) return fN;
 676   Int_t b=0, e=fN-1, m=(b+e)/2;
 677   for (; b<e; m=(b+e)/2) {
 678     if (z > fClusters[m]->GetZ()) b=m+1;
 679     else e=m;
 680   }
 681   return m;
 682 }
 683
 684 void AliITStrackerV2::AliITSlayer::
 685 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
 686   //--------------------------------------------------------------------
 687   // This function sets the "window"
 688   //--------------------------------------------------------------------
 689   fI=FindClusterIndex(zmin); fZmax=zmax;
 690   Double_t circle=2*TMath::Pi()*fR;
 691   if (ymax>circle) { ymax-=circle; ymin-=circle; }
 692   fYmin=ymin; fYmax=ymax;
 693 }
 694
 695 const AliITSclusterV2 *AliITStrackerV2::AliITSlayer::GetNextCluster(Int_t &ci){
 696   //--------------------------------------------------------------------
 697   // This function returns clusters within the "window"
 698   //--------------------------------------------------------------------
 699   const AliITSclusterV2 *cluster=0;
 700   for (Int_t i=fI; i<fN; i++) {
 701     const AliITSclusterV2 *c=fClusters[i];
 702     if (c->GetZ() > fZmax) break;
 703     if (c->IsUsed()) continue;
 704     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
 705     Double_t y=fR*det.GetPhi() + c->GetY();
 706
 707     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
 708     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
 709
 710     if (y<fYmin) continue;
 711     if (y>fYmax) continue;
 712     cluster=c; ci=i;
 713     fI=i+1;
 714     break;
 715   }
 716   return cluster;
 717 }
 718
 719 Int_t AliITStrackerV2::AliITSlayer::
 720 FindDetectorIndex(Double_t phi, Double_t z) const {
 721   //--------------------------------------------------------------------
 722   //This function finds the detector crossed by the track
 723   //--------------------------------------------------------------------
 724   Double_t dphi=phi-fPhiOffset;
 725   if      (dphi <  0) dphi += 2*TMath::Pi();
 726   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
 727   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
 728
 729   Double_t dz=fZOffset-z;
 730   Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
 731
 732   if (np>=fNladders) np-=fNladders;
 733   if (np<0)          np+=fNladders;
 734
 735 #ifdef DEBUG
 736 cout<<np<<' '<<nz<<endl;
 737 #endif
 738
 739   return np*fNdetectors + nz;
 740 }
 741
 742 Double_t
 743 AliITStrackerV2::AliITSlayer::GetThickness(Double_t phi, Double_t z) const {
 744   //--------------------------------------------------------------------
 745   //This function returns the thickness of this layer
 746   //--------------------------------------------------------------------
 747   //-pi<phi<+pi
 748   if (3 <fR&&fR<8 ) return 1.1*0.096;
 749   if (13<fR&&fR<26) return 1.1*0.088;
 750   if (37<fR&&fR<41) return 1.1*0.085;
 751   return 1.1*0.081;
 752 }
 753
 754
 755 Double_t AliITStrackerV2::GetEffectiveThickness(Double_t phi,Double_t z) const
 756 {
 757   //--------------------------------------------------------------------
 758   //Returns the thickness between the current layer and the vertex
 759   //--------------------------------------------------------------------
 760   //-pi<phi<+pi
 761   Double_t d=0.1;
 762
 763   Double_t xn=fLayers[fI].GetR();
 764   for (Int_t i=0; i<fI; i++) {
 765     Double_t xi=fLayers[i].GetR();
 766     d+=fLayers[i].GetThickness(phi,z)*xi*xi;
 767   }
 768
 769   if (fI>1) {
 770     Double_t xi=0.5*(fLayers[1].GetR()+fLayers[2].GetR());
 771     d+=0.034*xi*xi;
 772   }
 773
 774   if (fI>3) {
 775     Double_t xi=0.5*(fLayers[3].GetR()+fLayers[4].GetR());
 776     d+=0.039*xi*xi;
 777   }
 778   return d/(xn*xn);
 779 }
 780
 781
 782
 783 Int_t AliITStrackerV2::AliITSlayer::InRoad() const {
 784   //--------------------------------------------------------------------
 785   // This function returns number of clusters within the "window"
 786   //--------------------------------------------------------------------
 787   Int_t ncl=0;
 788   for (Int_t i=fI; i<fN; i++) {
 789     const AliITSclusterV2 *c=fClusters[i];
 790     if (c->GetZ() > fZmax) break;
 791     //if (c->IsUsed()) continue;
 792     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
 793     Double_t y=fR*det.GetPhi() + c->GetY();
 794
 795     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
 796     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
 797
 798     if (y<fYmin) continue;
 799     if (y>fYmax) continue;
 800     ncl++;
 801   }
 802   return ncl;
 803 }
 804
 805
 806