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