HLT/ITS/tracking/AliITStrackerHLT.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, 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 /* $Id: AliITStrackerHLT.cxx 32466 2009-05-20 07:51:56Z hristov $ */
  16
  17 //-------------------------------------------------------------------------
  18 //               Implementation of the ITS tracker class
  19 //    It reads AliITSRecPoint clusters and creates AliHLTITSTrack tracks
  20 //                   and fills with them the ESD
  21 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  22 //          Current support and development:
  23 //                     Andrea Dainese, andrea.dainese@lnl.infn.it
  24 //     dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  25 //     Params moved to AliITSRecoParam by: Andrea Dainese, INFN
  26 //     Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
  27 //-------------------------------------------------------------------------
  28
  29 //#include <TMatrixD.h>
  30 #include <TTree.h>
  31 #include <TDatabasePDG.h>
  32 #include <TString.h>
  33 #include <TRandom.h>
  34 #include <TTreeStream.h>
  35
  36
  37 #include "AliLog.h"
  38 #include "AliITSCalibrationSPD.h"
  39 #include "AliITSCalibrationSDD.h"
  40 #include "AliITSCalibrationSSD.h"
  41 #include "AliCDBEntry.h"
  42 #include "AliCDBManager.h"
  43 #include "AliAlignObj.h"
  44 #include "AliTrackPointArray.h"
  45 #include "AliESDVertex.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDtrack.h"
  48 #include "AliV0.h"
  49 #include "AliHelix.h"
  50 #include "AliITSChannelStatus.h"
  51 #include "AliITSDetTypeRec.h"
  52 #include "AliITSRecPoint.h"
  53 #include "AliITSgeomTGeo.h"
  54 #include "AliITSReconstructor.h"
  55 #include "AliITSClusterParam.h"
  56 #include "AliITSsegmentation.h"
  57 #include "AliITSCalibration.h"
  58 #include "AliITSV0Finder.h"
  59 #include "AliITStrackerHLT.h"
  60 //#include "AliHLTTPCCATrackParam.h"
  61 #include "AliHLTVertexer.h"
  62
  63
  64 ClassImp(AliITStrackerHLT)
  65
  66 Bool_t AliITStrackerHLT::TransportToX( AliExternalTrackParam *t, double x ) const
  67 {
  68   return t->PropagateTo( x, t->GetBz() );
  69 }
  70
  71 Bool_t AliITStrackerHLT::TransportToPhiX( AliExternalTrackParam *t, double phi, double x ) const
  72 {
  73   return t->Propagate( phi, x, t->GetBz() );
  74 }
  75
  76
  77 //------------------------------------------------------------------------
  78 Int_t AliITStrackerHLT::UpdateMI(AliHLTITSTrack* track, const AliITSRecPoint* cl,Double_t /*chi2*/,Int_t index) const
  79 {
  80   //
  81   // Update ITS track
  82   //
  83
  84   if (cl->GetQ()<=0) return 0;  // ingore the "virtual" clusters
  85
  86   Int_t layer = (index & 0xf0000000) >> 28;
  87
  88   // Take into account the mis-alignment (bring track to cluster plane)
  89
  90   Double_t xTrOrig=track->GetX();
  91   if (!TransportToX( track, xTrOrig + cl->GetX() ) ) return 0;
  92
  93   Double_t err2Y, err2Z;
  94
  95   GetClusterErrors2( layer, cl, track, err2Y, err2Z );
  96
  97   Double_t p[2]={ cl->GetY(), cl->GetZ()};
  98   Double_t cov[3]={err2Y, 0., err2Z};
  99
 100   Int_t updated = 1;
 101   //if( layer!=2 && layer!=3 )
 102   updated = track->AliExternalTrackParam::Update(p,cov);
 103
 104   int n = track->GetNumberOfClusters();
 105   track->SetClusterIndex(n,index);
 106   track->SetNumberOfClusters(n+1);
 107
 108   return updated;
 109 }
 110
 111
 112 AliHLTITSLayer AliITStrackerHLT::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
 113
 114 AliITStrackerHLT::AliITStrackerHLT():AliTracker(),
 115 fEsd(0),
 116 fUseTGeo(3),
 117 fxOverX0Pipe(-1.),
 118 fxTimesRhoPipe(-1.),
 119 fDebugStreamer(0),
 120 fITSChannelStatus(0),
 121 fkDetTypeRec(0)
 122 {
 123   //Default constructor
 124   Int_t i;
 125   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
 126   for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 127   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 128 }
 129 //------------------------------------------------------------------------
 130 AliITStrackerHLT::AliITStrackerHLT(const Char_t *geom) : AliTracker(),
 131 fEsd(0),
 132 fUseTGeo(3),
 133 fxOverX0Pipe(-1.),
 134 fxTimesRhoPipe(-1.),
 135 fDebugStreamer(0),
 136 fITSChannelStatus(0),
 137 fkDetTypeRec(0)
 138 {
 139   //--------------------------------------------------------------------
 140   //This is the AliITStrackerHLT constructor
 141   //--------------------------------------------------------------------
 142   if (geom) {
 143     AliWarning("\"geom\" is actually a dummy argument !");
 144   }
 145
 146   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 147     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
 148     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
 149
 150     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
 151     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 152     Double_t poff=TMath::ATan2(y,x);
 153     Double_t zoff=z;
 154     Double_t r=TMath::Sqrt(x*x + y*y);
 155
 156     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 157     r += TMath::Sqrt(x*x + y*y);
 158     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 159     r += TMath::Sqrt(x*x + y*y);
 160     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 161     r += TMath::Sqrt(x*x + y*y);
 162     r*=0.25;
 163
 164     new (fgLayers+i-1) AliHLTITSLayer(r,poff,zoff,nlad,ndet);
 165
 166     for (Int_t j=1; j<nlad+1; j++) {
 167       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 168         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 169         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 170         m.Multiply(tm);
 171         Double_t txyz[3]={0.};
 172         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 173         m.LocalToMaster(txyz,xyz);
 174         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 175         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 176
 177         if (phi<0) phi+=TMath::TwoPi();
 178         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 179
 180         AliHLTITSDetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 181         new(&det) AliHLTITSDetector(r,phi);
 182         // compute the real radius (with misalignment)
 183         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
 184         mmisal.Multiply(tm);
 185         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 186         mmisal.LocalToMaster(txyz,xyz);
 187         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 188         det.SetRmisal(rmisal);
 189
 190       } // end loop on detectors
 191     } // end loop on ladders
 192   } // end loop on layers
 193
 194
 195
 196   Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
 197                      AliITSReconstructor::GetRecoParam()->GetYVdef(),
 198                      AliITSReconstructor::GetRecoParam()->GetZVdef()};
 199   Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
 200                      AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
 201                      AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
 202   SetVertex(xyzVtx,ersVtx);
 203
 204   // store positions of centre of SPD modules (in z)
 205   Double_t tr[3];
 206   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 207   fSPDdetzcentre[0] = tr[2];
 208   AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 209   fSPDdetzcentre[1] = tr[2];
 210   AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 211   fSPDdetzcentre[2] = tr[2];
 212   AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 213   fSPDdetzcentre[3] = tr[2];
 214
 215   fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
 216   if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
 217     AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
 218     fUseTGeo = 3;
 219   }
 220
 221   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 222   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 223
 224   fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 225
 226 }
 227 //------------------------------------------------------------------------
 228 AliITStrackerHLT::AliITStrackerHLT(const AliITStrackerHLT &tracker):AliTracker(tracker),
 229 fEsd(tracker.fEsd),
 230 fUseTGeo(tracker.fUseTGeo),
 231 fxOverX0Pipe(tracker.fxOverX0Pipe),
 232 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
 233 fDebugStreamer(tracker.fDebugStreamer),
 234 fITSChannelStatus(tracker.fITSChannelStatus),
 235 fkDetTypeRec(tracker.fkDetTypeRec)
 236 {
 237   //Copy constructor
 238   Int_t i;
 239   for(i=0;i<4;i++) {
 240     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
 241   }
 242   for(i=0;i<6;i++) {
 243     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
 244     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
 245   }
 246   for(i=0;i<2;i++) {
 247     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
 248     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
 249   }
 250 }
 251 //------------------------------------------------------------------------
 252 AliITStrackerHLT & AliITStrackerHLT::operator=(const AliITStrackerHLT &tracker){
 253   //Assignment operator
 254   this->~AliITStrackerHLT();
 255   new(this) AliITStrackerHLT(tracker);
 256   return *this;
 257 }
 258 //------------------------------------------------------------------------
 259 AliITStrackerHLT::~AliITStrackerHLT()
 260 {
 261   //
 262   //destructor
 263   //
 264   if (fDebugStreamer) {
 265     //fDebugStreamer->Close();
 266     delete fDebugStreamer;
 267   }
 268   if(fITSChannelStatus) delete fITSChannelStatus;
 269 }
 270
 271
 272 //------------------------------------------------------------------------
 273 void AliITStrackerHLT::ReadBadFromDetTypeRec() {
 274   //--------------------------------------------------------------------
 275   //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
 276   //i.e. from OCDB
 277   //--------------------------------------------------------------------
 278
 279   if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
 280
 281   Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
 282
 283   if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
 284
 285   // ITS channels map
 286   if(fITSChannelStatus) delete fITSChannelStatus;
 287   fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
 288
 289   // ITS detectors and chips
 290   Int_t i=0,j=0,k=0,ndet=0;
 291   for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 292     Int_t nBadDetsPerLayer=0;
 293     ndet=AliITSgeomTGeo::GetNDetectors(i);
 294     for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
 295       for (k=1; k<ndet+1; k++) {
 296         AliHLTITSDetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 297         det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
 298         if(det.IsBad()) {nBadDetsPerLayer++;}
 299       } // end loop on detectors
 300     } // end loop on ladders
 301     Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
 302   } // end loop on layers
 303
 304   return;
 305 }
 306 //------------------------------------------------------------------------
 307 Int_t AliITStrackerHLT::LoadClusters(TTree *cTree) {
 308   //--------------------------------------------------------------------
 309   //This function loads ITS clusters
 310   //--------------------------------------------------------------------
 311   TBranch *branch=cTree->GetBranch("ITSRecPoints");
 312   if (!branch) {
 313     Error("LoadClusters"," can't get the branch !\n");
 314     return 1;
 315   }
 316
 317   static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 318   branch->SetAddress(&clusters);
 319
 320   Int_t i=0,j=0,ndet=0;
 321   Int_t detector=0;
 322   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 323     ndet=fgLayers[i].GetNdetectors();
 324     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 325     for (; j<jmax; j++) {
 326       if (!cTree->GetEvent(j)) continue;
 327       Int_t ncl=clusters->GetEntriesFast();
 328       SignDeltas(clusters,GetZ());
 329
 330       while (ncl--) {
 331         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 332         detector=c->GetDetectorIndex();
 333
 334         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 335
 336         fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 337       }
 338       clusters->Delete();
 339       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 340       // zwindow cm from the dead zone
 341       if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
 342         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
 343           Int_t lab[4]   = {0,0,0,detector};
 344           Int_t info[3]  = {0,0,i};
 345           Float_t q      = 0.; // this identifies virtual clusters
 346           Float_t hit[5] = {xdead,
 347                             0.,
 348                             AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
 349                             AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
 350                             q};
 351           Bool_t local   = kTRUE;
 352           Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
 353           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 354           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 355             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 356           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 357           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 358             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 359           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 360           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 361             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 362           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 363           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 364             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 365           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 366           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 367             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 368           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
 369           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 370             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 371         }
 372       } // "virtual" clusters in SPD
 373
 374     }
 375     //
 376     fgLayers[i].ResetRoad(); //road defined by the cluster density
 377     fgLayers[i].SortClusters();
 378   }
 379
 380   dummy.Clear();
 381
 382   return 0;
 383 }
 384 //------------------------------------------------------------------------
 385 void AliITStrackerHLT::UnloadClusters() {
 386   //--------------------------------------------------------------------
 387   //This function unloads ITS clusters
 388   //--------------------------------------------------------------------
 389   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
 390 }
 391
 392
 393
 394 //------------------------------------------------------------------------
 395 Int_t AliITStrackerHLT::CorrectForTPCtoITSDeadZoneMaterial(AliHLTITSTrack *t) {
 396   //--------------------------------------------------------------------
 397   // Correction for the material between the TPC and the ITS
 398   //--------------------------------------------------------------------
 399   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction
 400       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
 401       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
 402       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
 403   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
 404       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
 405       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
 406       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
 407   } else {
 408     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
 409     return 0;
 410   }
 411
 412   return 1;
 413 }
 414
 415 #include "TStopwatch.h"
 416
 417 //------------------------------------------------------------------------
 418 Int_t AliITStrackerHLT::Clusters2Tracks(AliESDEvent *event) {
 419   //--------------------------------------------------------------------
 420   // This functions reconstructs ITS tracks
 421   // The clusters must be already loaded !
 422   //--------------------------------------------------------------------
 423   std::cout<<"\n\n ITS starts...\n"<<std::endl;
 424   TStopwatch timer;
 425
 426   fEsd = event;         // store pointer to the esd
 427
 428   {/* Read ESD tracks */
 429
 430     Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 431
 432     for( int itr=0; itr<event->GetNumberOfTracks(); itr++ ){
 433
 434       AliESDtrack *esdTrack = event->GetTrack(itr);
 435
 436       if ((esdTrack->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 437       if (esdTrack->GetStatus()&AliESDtrack::kTPCout) continue;
 438       if (esdTrack->GetStatus()&AliESDtrack::kITSin) continue;
 439       if (esdTrack->GetKinkIndex(0)>0) continue;   //kink daughter
 440       AliHLTITSTrack tMI(*esdTrack);
 441       AliHLTITSTrack *t = &tMI;
 442
 443       if (tMI.GetMass()<0.9*pimass) t->SetMass(pimass);
 444       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 445
 446       if (!CorrectForTPCtoITSDeadZoneMaterial(t))  continue;
 447
 448       Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 449
 450       FollowProlongationTree(t);
 451       int nclu=0;
 452       for(Int_t i=0; i<6; i++) {
 453         if( t->GetClusterIndex(i)>=0 ) nclu++;
 454       }
 455       cout<<"N assigned ITS clusters = "<<nclu<<std::endl;
 456       if( nclu>0 ){
 457         t->SetLabel(-1);//tpcLabel);
 458         t->SetFakeRatio(1.);
 459         CookLabel(t,0.); //For comparison only
 460         cout<<"label = "<<t->GetLabel()<<" / "<<tpcLabel<<endl;
 461         TransportToX(t, 0 );
 462         cout<<"\n fill track : parameters at "<<t->GetX()<<": "<< TMath::Sqrt(TMath::Abs(t->GetSigmaY2()))<<" "<< TMath::Sqrt(TMath::Abs(t->GetSigmaY2()))<<endl;
 463         //t->Print();
 464         UpdateESDtrack(t,AliESDtrack::kITSin);
 465       }
 466     }
 467   } /* End Read ESD tracks */
 468
 469   AliHLTVertexer vertexer;
 470   vertexer.SetESD( event );
 471   vertexer.FindPrimaryVertex();
 472   vertexer.FindV0s();
 473
 474
 475   timer.Stop();
 476   static double totalTime = 0;
 477   static int nEvnts = 0;
 478   totalTime+=timer.CpuTime();
 479   nEvnts++;
 480   std::cout<<"\n\n ITS tracker time = "<<totalTime/nEvnts<<" [s/ev]  for "<<nEvnts<<" events\n\n "<<std::endl;
 481   return 0;
 482 }
 483
 484
 485 //------------------------------------------------------------------------
 486 Int_t AliITStrackerHLT::PropagateBack(AliESDEvent * /*event*/) {
 487   return 0;
 488 }
 489
 490 //------------------------------------------------------------------------
 491 Int_t AliITStrackerHLT::RefitInward(AliESDEvent * /*event*/ ) {
 492   return 0;
 493 }
 494 //------------------------------------------------------------------------
 495 AliCluster *AliITStrackerHLT::GetCluster(Int_t index) const {
 496   //--------------------------------------------------------------------
 497   //       Return pointer to a given cluster
 498   //--------------------------------------------------------------------
 499   Int_t l=(index & 0xf0000000) >> 28;
 500   Int_t c=(index & 0x0fffffff) >> 00;
 501   return fgLayers[l].GetCluster(c);
 502 }
 503 //------------------------------------------------------------------------
 504 Bool_t AliITStrackerHLT::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 505   //--------------------------------------------------------------------
 506   // Get track space point with index i
 507   //--------------------------------------------------------------------
 508
 509   Int_t l=(index & 0xf0000000) >> 28;
 510   Int_t c=(index & 0x0fffffff) >> 00;
 511   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 512   Int_t idet = cl->GetDetectorIndex();
 513
 514   Float_t xyz[3];
 515   Float_t cov[6];
 516   cl->GetGlobalXYZ(xyz);
 517   cl->GetGlobalCov(cov);
 518   p.SetXYZ(xyz, cov);
 519   p.SetCharge(cl->GetQ());
 520   p.SetDriftTime(cl->GetDriftTime());
 521   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 522   switch (l) {
 523   case 0:
 524     iLayer = AliGeomManager::kSPD1;
 525     break;
 526   case 1:
 527     iLayer = AliGeomManager::kSPD2;
 528     break;
 529   case 2:
 530     iLayer = AliGeomManager::kSDD1;
 531     break;
 532   case 3:
 533     iLayer = AliGeomManager::kSDD2;
 534     break;
 535   case 4:
 536     iLayer = AliGeomManager::kSSD1;
 537     break;
 538   case 5:
 539     iLayer = AliGeomManager::kSSD2;
 540     break;
 541   default:
 542     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 543     break;
 544   };
 545   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 546   p.SetVolumeID((UShort_t)volid);
 547   return kTRUE;
 548 }
 549 //------------------------------------------------------------------------
 550 Bool_t AliITStrackerHLT::GetTrackPointTrackingError(Int_t index,
 551                         AliTrackPoint& p, const AliESDtrack *t) {
 552   //--------------------------------------------------------------------
 553   // Get track space point with index i
 554   // (assign error estimated during the tracking)
 555   //--------------------------------------------------------------------
 556
 557   Int_t l=(index & 0xf0000000) >> 28;
 558   Int_t c=(index & 0x0fffffff) >> 00;
 559   const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 560   Int_t idet = cl->GetDetectorIndex();
 561
 562   const AliHLTITSDetector &det=fgLayers[l].GetDetector(idet);
 563
 564   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
 565   Float_t detxy[2];
 566   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
 567   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
 568   Double_t alpha = t->GetAlpha();
 569   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
 570   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
 571   phi += alpha-det.GetPhi();
 572   Float_t tgphi = TMath::Tan(phi);
 573
 574   Float_t tgl = t->GetTgl(); // tgl about const along track
 575   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
 576
 577   Float_t errlocalx,errlocalz;
 578   Bool_t addMisalErr=kFALSE;
 579   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errlocalx,errlocalz,addMisalErr);
 580
 581   Float_t xyz[3];
 582   Float_t cov[6];
 583   cl->GetGlobalXYZ(xyz);
 584   //  cl->GetGlobalCov(cov);
 585   Float_t pos[3] = {0.,0.,0.};
 586   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errlocalx*errlocalx,errlocalz*errlocalz,0);
 587   tmpcl.GetGlobalCov(cov);
 588
 589   p.SetXYZ(xyz, cov);
 590   p.SetCharge(cl->GetQ());
 591   p.SetDriftTime(cl->GetDriftTime());
 592
 593   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 594   switch (l) {
 595   case 0:
 596     iLayer = AliGeomManager::kSPD1;
 597     break;
 598   case 1:
 599     iLayer = AliGeomManager::kSPD2;
 600     break;
 601   case 2:
 602     iLayer = AliGeomManager::kSDD1;
 603     break;
 604   case 3:
 605     iLayer = AliGeomManager::kSDD2;
 606     break;
 607   case 4:
 608     iLayer = AliGeomManager::kSSD1;
 609     break;
 610   case 5:
 611     iLayer = AliGeomManager::kSSD2;
 612     break;
 613   default:
 614     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 615     break;
 616   };
 617   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 618
 619   p.SetVolumeID((UShort_t)volid);
 620   return kTRUE;
 621 }
 622
 623
 624 //------------------------------------------------------------------------
 625 void AliITStrackerHLT::FollowProlongationTree(AliHLTITSTrack * track )
 626 {
 627     //cout<<endl;
 628   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
 629     //cout<<"\nLayer "<<ilayer<<endl;
 630
 631     AliHLTITSLayer &layer=fgLayers[ilayer];
 632
 633     //cout<<" shield material.. "<<endl;
 634
 635     // material between SSD and SDD, SDD and SPD
 636     if (ilayer==3 && !CorrectForShieldMaterial(track,"SDD","inward")) continue;
 637     if (ilayer==1 && !CorrectForShieldMaterial(track,"SPD","inward")) continue;
 638
 639     int idet;
 640
 641     {
 642       // propagate to the layer radius
 643
 644       Double_t r = layer.GetR(), phi,z;
 645       Double_t xToGo;
 646       //cout<<" propagate to layer R= "<<r<<" .."<<endl;
 647       if( !track->GetLocalXat(r,xToGo) ) continue;
 648       if( !TransportToX(track, xToGo) ) continue;
 649
 650       // detector number
 651
 652       if (!track->GetPhiZat(r,phi,z)) continue;
 653       idet=layer.FindDetectorIndex(phi,z);
 654       //cout<<" detector number = "<<idet<<endl;
 655    }
 656
 657
 658     //cout<<" correct for the layer material .. "<<endl;
 659
 660     // correct for the layer material
 661     {
 662       Double_t trackGlobXYZ1[3];
 663       if (!track->GetXYZ(trackGlobXYZ1)) continue;
 664       CorrectForLayerMaterial(track,ilayer,trackGlobXYZ1,"inward");
 665     }
 666
 667     // track outside layer acceptance in z
 668
 669     if( idet<0 ) continue;
 670
 671     // propagate to the intersection with the detector plane
 672
 673     //cout<<" propagate to the intersection with the detector .. "<<endl;
 674
 675     const AliHLTITSDetector &det=layer.GetDetector( idet );
 676     if (!TransportToPhiX( track, det.GetPhi(), det.GetR() ) ) continue;
 677
 678     // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
 679
 680     // road in global (rphi,z) [i.e. in tracking ref. system]
 681
 682     Double_t zmin,zmax,ymin,ymax;
 683
 684     //cout<<" ComputeRoad .. "<<endl;
 685     if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
 686
 687     //cout<<" Road: y["<<ymin<<","<<ymax<<"], z["<<zmin<<","<<zmax<<"] "<<endl;
 688
 689     // select clusters in road
 690
 691     //cout<<" SelectClusters .. "<<endl;
 692     layer.SelectClusters(zmin,zmax,ymin,ymax);
 693
 694     // Define criteria for track-cluster association
 695
 696     Double_t msz = track->GetSigmaZ2() +
 697       AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
 698       AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
 699       AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
 700
 701     Double_t msy = track->GetSigmaY2() +
 702       AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
 703       AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
 704       AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
 705
 706      msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
 707      msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
 708
 709      msz = 1./msz; // 1/RoadZ^2
 710      msy = 1./msy; // 1/RoadY^2
 711
 712      //
 713      // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
 714      //
 715
 716      const AliITSRecPoint *cl=0;
 717      Int_t clidx=-1;
 718      Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
 719      Bool_t deadzoneSPD=kFALSE;
 720
 721      // check if the road contains a dead zone
 722      Bool_t noClusters = !layer.GetNextCluster(clidx,kTRUE);
 723
 724      Double_t dz=0.5*(zmax-zmin);
 725      Double_t dy=0.5*(ymax-ymin);
 726
 727      Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,noClusters);
 728
 729      // create a prolongation without clusters (check also if there are no clusters in the road)
 730
 731      if (dead==1) { // dead zone at z=0,+-7cm in SPD
 732        deadzoneSPD=kTRUE;
 733      }
 734
 735      clidx=-1;
 736
 737      //cout<<" loop over clusters in the road .. "<<endl;
 738
 739      // loop over clusters in the road
 740      const AliITSRecPoint *bestCluster=0;
 741      double bestChi2 = 1.e10;
 742      AliHLTITSTrack bestTrack( *track );
 743      int bestIdx = -1;
 744      while ((cl=layer.GetNextCluster(clidx))!=0) {
 745        //cout<<" cluster: "<<cl->GetX()<<" "<<cl->GetY()<<" "<<cl->GetZ()<<endl;
 746        AliHLTITSTrack t(*track);
 747        if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
 748
 749        Int_t idetc=cl->GetDetectorIndex();
 750
 751        //cout<<" cluster detector: "<<idetc<<endl;
 752
 753        if ( idet !=idetc ) { // new cluster's detector
 754          const AliHLTITSDetector &detc=layer.GetDetector(idetc);
 755          if (!TransportToPhiX( track, detc.GetPhi(),detc.GetR()) ) continue;
 756          t = *track;
 757          idet = idetc;
 758        }
 759
 760        // take into account misalignment (bring track to real detector plane)
 761
 762        if (!TransportToX( &t, t.GetX() + cl->GetX() ) ) continue;
 763        double chi2 = ( (t.GetZ()-cl->GetZ())*(t.GetZ()-cl->GetZ())*msz +
 764                        (t.GetY()-cl->GetY())*(t.GetY()-cl->GetY())*msy   );
 765        //cout<<" chi2="<<chi2<<endl;
 766        if ( chi2 < bestChi2 ){
 767          bestChi2 = chi2;
 768          bestCluster = cl;
 769          bestTrack = t;
 770          bestIdx = clidx;
 771          continue;
 772        }
 773      }
 774
 775      //cout<<" best chi2= "<<bestChi2<<endl;
 776
 777      if( bestCluster && bestChi2 <=1. ){
 778
 779        //cout<<" cluster found "<<endl;
 780        *track = bestTrack;
 781
 782
 783       // calculate track-clusters chi2
 784        chi2trkcl = GetPredictedChi2MI(track,bestCluster,ilayer);
 785        //cout<<" track-clusters chi2 = "<<chi2trkcl<<endl;
 786        //cout<<" max chi2 = "<<AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)<<endl;
 787
 788        // chi2 cut
 789        if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
 790          if (bestCluster->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
 791          //cout<<"set index.."<<endl;
 792          //cout<<"set index ok"<<endl;
 793          if (bestCluster->GetQ()!=0) { // real cluster
 794            //cout<<" UpdateMI ... "<<endl;
 795            if (!UpdateMI(track,bestCluster,chi2trkcl,(ilayer<<28)+bestIdx)) {
 796              continue;
 797            }
 798          }
 799        }
 800      }
 801      //cout<<" goto next layer "<<endl;
 802
 803   }
 804 }
 805
 806
 807 //------------------------------------------------------------------------
 808 AliHLTITSLayer & AliITStrackerHLT::GetLayer(Int_t layer) const
 809 {
 810   //--------------------------------------------------------------------
 811   //
 812   //
 813   return fgLayers[layer];
 814 }
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828 //------------------------------------------------------------------------
 829 void AliITStrackerHLT::CookLabel(AliHLTITSTrack *track,Float_t wrong) const {
 830   //--------------------------------------------------------------------
 831   //This function "cooks" a track label. If label<0, this track is fake.
 832   //--------------------------------------------------------------------
 833   Int_t tpcLabel=-1;
 834
 835   if ( track->GetESDtrack())   tpcLabel =  TMath::Abs(track->GetESDtrack()->GetTPCLabel());
 836
 837    Int_t nwrong=0;
 838    cout<<"cook label: nclu = "<<track->GetNumberOfClusters()<<endl;
 839    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
 840      Int_t cindex = track->GetClusterIndex(i);
 841      //Int_t l=(cindex & 0xf0000000) >> 28;
 842      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
 843      Int_t isWrong=1;
 844      for (Int_t ind=0;ind<3;ind++){
 845        if (tpcLabel>0)
 846          if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
 847        AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
 848      }
 849       nwrong+=isWrong;
 850    }
 851    Int_t nclusters = track->GetNumberOfClusters();
 852    if (nclusters > 0) //PH Some tracks don't have any cluster
 853      track->SetFakeRatio(double(nwrong)/double(nclusters));
 854    cout<<"fake ratio = "<<track->GetFakeRatio()<<endl;
 855    if (tpcLabel>0){
 856      if (track->GetFakeRatio()>wrong) track->SetLabel(-tpcLabel);
 857      else
 858        track->SetLabel(tpcLabel);
 859    }
 860    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
 861
 862 }
 863
 864
 865
 866 void AliITStrackerHLT::GetClusterErrors2( Int_t layer, const AliITSRecPoint *cluster, AliHLTITSTrack* track, double &err2Y, double &err2Z ) const
 867 {
 868   Float_t erry,errz;
 869   Float_t theta = track->GetTgl();
 870   Float_t phi   = track->GetSnp();
 871   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
 872   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz);
 873   err2Y = erry*erry;
 874   err2Z = errz*errz;
 875 }
 876
 877
 878 //------------------------------------------------------------------------
 879 Double_t AliITStrackerHLT::GetPredictedChi2MI(AliHLTITSTrack* track, const AliITSRecPoint *cluster,Int_t layer)
 880 {
 881   //
 882   // Compute predicted chi2
 883   //
 884
 885   AliHLTITSTrack t(*track);
 886   if (!t.Propagate( t.GetX()+cluster->GetX())) return 1000.;
 887
 888   Double_t  err2Y,err2Z;
 889   GetClusterErrors2( layer, cluster, &t, err2Y, err2Z );
 890
 891   Double_t chi2 = t.GetPredictedChi2(cluster->GetY(),cluster->GetZ(),err2Y,err2Z);
 892
 893   Float_t ny,nz;
 894   Float_t theta = track->GetTgl();
 895   Float_t phi   = track->GetSnp();
 896   phi = TMath::Abs(phi)*TMath::Sqrt(1./((1.-phi)*(1.+phi)));
 897   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
 898   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
 899   if (delta>1){
 900     chi2+=0.5*TMath::Min(delta/2,2.);
 901     chi2+=2.*cluster->GetDeltaProbability();
 902   }
 903   return chi2;
 904 }
 905
 906
 907 //------------------------------------------------------------------------
 908 void AliITStrackerHLT::SignDeltas(const TObjArray *clusterArray, Float_t vz)
 909 {
 910   //
 911   // Clusters from delta electrons?
 912   //
 913   Int_t entries = clusterArray->GetEntriesFast();
 914   if (entries<4) return;
 915   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
 916   Int_t layer = cluster->GetLayer();
 917   if (layer>1) return;
 918   Int_t index[10000];
 919   Int_t ncandidates=0;
 920   Float_t r = (layer>0)? 7:4;
 921   //
 922   for (Int_t i=0;i<entries;i++){
 923     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
 924     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
 925     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
 926     index[ncandidates] = i;  //candidate to belong to delta electron track
 927     ncandidates++;
 928     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
 929       cl0->SetDeltaProbability(1);
 930     }
 931   }
 932   //
 933   //
 934   //
 935   for (Int_t i=0;i<ncandidates;i++){
 936     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
 937     if (cl0->GetDeltaProbability()>0.8) continue;
 938     //
 939     Int_t ncl = 0;
 940     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
 941     sumy=sumz=sumy2=sumyz=sumw=0.0;
 942     for (Int_t j=0;j<ncandidates;j++){
 943       if (i==j) continue;
 944       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
 945       //
 946       Float_t dz = cl0->GetZ()-cl1->GetZ();
 947       Float_t dy = cl0->GetY()-cl1->GetY();
 948       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
 949         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
 950         y[ncl] = cl1->GetY();
 951         z[ncl] = cl1->GetZ();
 952         sumy+= y[ncl]*weight;
 953         sumz+= z[ncl]*weight;
 954         sumy2+=y[ncl]*y[ncl]*weight;
 955         sumyz+=y[ncl]*z[ncl]*weight;
 956         sumw+=weight;
 957         ncl++;
 958       }
 959     }
 960     if (ncl<4) continue;
 961     Float_t det = sumw*sumy2  - sumy*sumy;
 962     Float_t delta=1000;
 963     if (TMath::Abs(det)>0.01){
 964       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
 965       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
 966       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
 967     }
 968     else{
 969       Float_t z0  = sumyz/sumy;
 970       delta = TMath::Abs(cl0->GetZ()-z0);
 971     }
 972     if ( delta<0.05) {
 973       cl0->SetDeltaProbability(1-20.*delta);
 974     }
 975   }
 976 }
 977 //------------------------------------------------------------------------
 978 void AliITStrackerHLT::UpdateESDtrack(AliHLTITSTrack* track, ULong_t flags) const
 979 {
 980   //
 981   // Update ESD track
 982   //
 983   track->UpdateESDtrack(flags);
 984   AliHLTITSTrack * oldtrack = (AliHLTITSTrack*)(track->GetESDtrack()->GetITStrack());
 985   if (oldtrack) delete oldtrack;
 986   track->GetESDtrack()->SetITStrack(new AliHLTITSTrack(*track));
 987 }
 988
 989
 990
 991
 992 //------------------------------------------------------------------------
 993 void AliITStrackerHLT::BuildMaterialLUT(TString material) {
 994   //--------------------------------------------------------------------
 995   // Fill a look-up table with mean material
 996   //--------------------------------------------------------------------
 997
 998   Int_t n=1000;
 999   Double_t mparam[7];
1000   Double_t point1[3],point2[3];
1001   Double_t phi,cosphi,sinphi,z;
1002   // 0-5 layers, 6 pipe, 7-8 shields
1003   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
1004   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
1005
1006   Int_t ifirst=0,ilast=0;
1007   if(material.Contains("Pipe")) {
1008     ifirst=6; ilast=6;
1009   } else if(material.Contains("Shields")) {
1010     ifirst=7; ilast=8;
1011   } else if(material.Contains("Layers")) {
1012     ifirst=0; ilast=5;
1013   } else {
1014     Error("BuildMaterialLUT","Wrong layer name\n");
1015   }
1016
1017   for(Int_t imat=ifirst; imat<=ilast; imat++) {
1018     Double_t param[5]={0.,0.,0.,0.,0.};
1019     for (Int_t i=0; i<n; i++) {
1020       phi = 2.*TMath::Pi()*gRandom->Rndm();
1021       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
1022       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
1023       point1[0] = rmin[imat]*cosphi;
1024       point1[1] = rmin[imat]*sinphi;
1025       point1[2] = z;
1026       point2[0] = rmax[imat]*cosphi;
1027       point2[1] = rmax[imat]*sinphi;
1028       point2[2] = z;
1029       AliTracker::MeanMaterialBudget(point1,point2,mparam);
1030       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
1031     }
1032     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
1033     if(imat<=5) {
1034       fxOverX0Layer[imat] = param[1];
1035       fxTimesRhoLayer[imat] = param[0]*param[4];
1036     } else if(imat==6) {
1037       fxOverX0Pipe = param[1];
1038       fxTimesRhoPipe = param[0]*param[4];
1039     } else if(imat==7) {
1040       fxOverX0Shield[0] = param[1];
1041       fxTimesRhoShield[0] = param[0]*param[4];
1042     } else if(imat==8) {
1043       fxOverX0Shield[1] = param[1];
1044       fxTimesRhoShield[1] = param[0]*param[4];
1045     }
1046   }
1047   /*
1048   printf("%s\n",material.Data());
1049   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
1050   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
1051   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
1052   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
1053   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
1054   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
1055   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
1056   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
1057   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
1058   */
1059   return;
1060 }
1061
1062 //------------------------------------------------------------------------
1063 Int_t AliITStrackerHLT::CorrectForPipeMaterial(AliHLTITSTrack *t,
1064                                               TString direction) {
1065   //-------------------------------------------------------------------
1066   // Propagate beyond beam pipe and correct for material
1067   // (material budget in different ways according to fUseTGeo value)
1068   // Add time if going outward (PropagateTo or PropagateToTGeo)
1069   //-------------------------------------------------------------------
1070
1071   // Define budget mode:
1072   // 0: material from AliITSRecoParam (hard coded)
1073   // 1: material from TGeo in one step (on the fly)
1074   // 2: material from lut
1075   // 3: material from TGeo in one step (same for all hypotheses)
1076   Int_t mode;
1077   switch(fUseTGeo) {
1078   case 0:
1079     mode=0;
1080     break;
1081   case 1:
1082     mode=1;
1083     break;
1084   case 2:
1085     mode=2;
1086     break;
1087   case 3:
1088     mode=3;
1089     break;
1090   case 4:
1091     mode=3;
1092     break;
1093   default:
1094     mode=0;
1095     break;
1096   }
1097
1098   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1099   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
1100   Double_t xToGo;
1101   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1102
1103   Double_t xOverX0,x0,lengthTimesMeanDensity;
1104
1105   switch(mode) {
1106   case 0:
1107     xOverX0 = AliITSRecoParam::GetdPipe();
1108     x0 = AliITSRecoParam::GetX0Be();
1109     lengthTimesMeanDensity = xOverX0*x0;
1110     lengthTimesMeanDensity *= dir;
1111     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1112     break;
1113   case 1:
1114     if (!t->PropagateToTGeo(xToGo,1)) return 0;
1115     break;
1116   case 2:
1117     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
1118     xOverX0 = fxOverX0Pipe;
1119     lengthTimesMeanDensity = fxTimesRhoPipe;
1120     lengthTimesMeanDensity *= dir;
1121     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1122     break;
1123   case 3:
1124     double xOverX0PipeTrks, xTimesRhoPipeTrks;
1125     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1126     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1127                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
1128     xOverX0PipeTrks = TMath::Abs(xOverX0)/angle;
1129     xTimesRhoPipeTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1130     xOverX0 = xOverX0PipeTrks;
1131     lengthTimesMeanDensity = xTimesRhoPipeTrks;
1132     lengthTimesMeanDensity *= dir;
1133     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1134     break;
1135   }
1136
1137   return 1;
1138 }
1139 //------------------------------------------------------------------------
1140 Int_t AliITStrackerHLT::CorrectForShieldMaterial(AliHLTITSTrack *t,
1141                                                 TString shield,
1142                                                 TString direction) {
1143   //-------------------------------------------------------------------
1144   // Propagate beyond SPD or SDD shield and correct for material
1145   // (material budget in different ways according to fUseTGeo value)
1146   // Add time if going outward (PropagateTo or PropagateToTGeo)
1147   //-------------------------------------------------------------------
1148
1149   // Define budget mode:
1150   // 0: material from AliITSRecoParam (hard coded)
1151   // 1: material from TGeo in steps of X cm (on the fly)
1152   //    X = AliITSRecoParam::GetStepSizeTGeo()
1153   // 2: material from lut
1154   // 3: material from TGeo in one step (same for all hypotheses)
1155   Int_t mode;
1156   switch(fUseTGeo) {
1157   case 0:
1158     mode=0;
1159     break;
1160   case 1:
1161     mode=1;
1162     break;
1163   case 2:
1164     mode=2;
1165     break;
1166   case 3:
1167     mode=3;
1168     break;
1169   case 4:
1170     mode=3;
1171     break;
1172   default:
1173     mode=0;
1174     break;
1175   }
1176
1177
1178   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1179   Double_t rToGo;
1180   Int_t    shieldindex=0;
1181   if (shield.Contains("SDD")) { // SDDouter
1182     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
1183     shieldindex=1;
1184   } else if (shield.Contains("SPD")) {        // SPDouter
1185     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
1186     shieldindex=0;
1187   } else {
1188     Error("CorrectForShieldMaterial"," Wrong shield name\n");
1189     return 0;
1190   }
1191   Double_t xToGo;
1192   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1193
1194   Double_t xOverX0,x0,lengthTimesMeanDensity;
1195   Int_t nsteps=1;
1196
1197   switch(mode) {
1198   case 0:
1199     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
1200     x0 = AliITSRecoParam::GetX0shield(shieldindex);
1201     lengthTimesMeanDensity = xOverX0*x0;
1202     lengthTimesMeanDensity *= dir;
1203     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1204     break;
1205   case 1:
1206     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
1207     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1208     break;
1209   case 2:
1210     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
1211     xOverX0 = fxOverX0Shield[shieldindex];
1212     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
1213     lengthTimesMeanDensity *= dir;
1214     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1215     break;
1216   case 3:
1217     if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
1218     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1219                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
1220     double xOverX0ShieldTrks = TMath::Abs(xOverX0)/angle;
1221     double xTimesRhoShieldTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1222     xOverX0 = xOverX0ShieldTrks;
1223     lengthTimesMeanDensity = xTimesRhoShieldTrks;
1224     lengthTimesMeanDensity *= dir;
1225     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1226     break;
1227   }
1228
1229   return 1;
1230 }
1231 //------------------------------------------------------------------------
1232 Int_t AliITStrackerHLT::CorrectForLayerMaterial(AliHLTITSTrack *t,
1233                                                Int_t layerindex,
1234                                                Double_t oldGlobXYZ[3],
1235                                                TString direction) {
1236   //-------------------------------------------------------------------
1237   // Propagate beyond layer and correct for material
1238   // (material budget in different ways according to fUseTGeo value)
1239   // Add time if going outward (PropagateTo or PropagateToTGeo)
1240   //-------------------------------------------------------------------
1241
1242   // Define budget mode:
1243   // 0: material from AliITSRecoParam (hard coded)
1244   // 1: material from TGeo in stepsof X cm (on the fly)
1245   //    X = AliITSRecoParam::GetStepSizeTGeo()
1246   // 2: material from lut
1247   // 3: material from TGeo in one step (same for all hypotheses)
1248   Int_t mode;
1249   switch(fUseTGeo) {
1250   case 0:
1251     mode=0;
1252     break;
1253   case 1:
1254     mode=1;
1255     break;
1256   case 2:
1257     mode=2;
1258     break;
1259   case 3:
1260     mode=3;
1261     break;
1262   case 4:
1263     mode=3;
1264     break;
1265   default:
1266     mode=0;
1267     break;
1268   }
1269
1270
1271   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
1272
1273   Double_t r=fgLayers[layerindex].GetR();
1274   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
1275
1276   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
1277   Double_t xToGo;
1278   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
1279
1280
1281   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
1282   Int_t nsteps=1;
1283
1284   // back before material (no correction)
1285   Double_t rOld,xOld;
1286   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
1287   if (!t->GetLocalXat(rOld,xOld)) return 0;
1288   if (!TransportToX( t, xOld)) return 0;
1289
1290   switch(mode) {
1291   case 0:
1292     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
1293     lengthTimesMeanDensity = xOverX0*x0;
1294     lengthTimesMeanDensity *= dir;
1295     // Bring the track beyond the material
1296     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1297     break;
1298   case 1:
1299     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
1300     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
1301     break;
1302   case 2:
1303     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
1304     xOverX0 = fxOverX0Layer[layerindex];
1305     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
1306     lengthTimesMeanDensity *= dir;
1307     // Bring the track beyond the material
1308     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1309     break;
1310   case 3:
1311     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
1312     if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
1313     Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
1314                                ((1.-t->GetSnp())*(1.+t->GetSnp())));
1315     double xOverX0LayerTrks = TMath::Abs(xOverX0)/angle;
1316     double xTimesRhoLayerTrks = TMath::Abs(lengthTimesMeanDensity)/angle;
1317     xOverX0 = xOverX0LayerTrks;
1318     lengthTimesMeanDensity = xTimesRhoLayerTrks;
1319     lengthTimesMeanDensity *= dir;
1320     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
1321     break;
1322   }
1323
1324
1325   return 1;
1326 }
1327
1328
1329 //------------------------------------------------------------------------
1330 Int_t AliITStrackerHLT::CheckDeadZone(AliHLTITSTrack *track,
1331                                      Int_t ilayer,Int_t idet,
1332                                      Double_t dz,Double_t dy,
1333                                      Bool_t noClusters) const {
1334   //-----------------------------------------------------------------
1335   // This method is used to decide whether to allow a prolongation
1336   // without clusters, because there is a dead zone in the road.
1337   // In this case the return value is > 0:
1338   // return 1: dead zone at z=0,+-7cm in SPD
1339   // return 2: all road is "bad" (dead or noisy) from the OCDB
1340   // return 3: something "bad" (dead or noisy) from the OCDB
1341   //-----------------------------------------------------------------
1342
1343   // check dead zones at z=0,+-7cm in the SPD
1344   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
1345     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
1346                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
1347                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
1348     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
1349                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
1350                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
1351     for (Int_t i=0; i<3; i++)
1352       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
1353         AliDebug(2,Form("crack SPD %d",ilayer));
1354         return 1;
1355       }
1356   }
1357
1358   // check bad zones from OCDB
1359   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
1360
1361   if (idet<0) return 0;
1362
1363   AliHLTITSDetector &det=fgLayers[ilayer].GetDetector(idet);
1364
1365   Int_t detType=-1;
1366   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
1367   if (ilayer==0 || ilayer==1) {        // ----------  SPD
1368     detType = 0;
1369   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
1370     detType = 1;
1371     detSizeFactorX *= 2.;
1372   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
1373     detType = 2;
1374   }
1375   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
1376   if (detType==2) segm->SetLayer(ilayer+1);
1377   Float_t detSizeX = detSizeFactorX*segm->Dx();
1378   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
1379
1380   // check if the road overlaps with bad chips
1381   Float_t xloc,zloc;
1382   LocalModuleCoord(ilayer,idet,track,xloc,zloc);
1383   Float_t zlocmin = zloc-dz;
1384   Float_t zlocmax = zloc+dz;
1385   Float_t xlocmin = xloc-dy;
1386   Float_t xlocmax = xloc+dy;
1387   Int_t chipsInRoad[100];
1388
1389   // check if road goes out of detector
1390   Bool_t touchNeighbourDet=kFALSE;
1391   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
1392   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
1393   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
1394   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
1395   AliDebug(2,Form("layer %d det %d zmim zmax %f %f xmin xmax %f %f   %f %f",ilayer,idet,zlocmin,zlocmax,xlocmin,xlocmax,detSizeZ,detSizeX));
1396
1397   // check if this detector is bad
1398   if (det.IsBad()) {
1399     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
1400     if(!touchNeighbourDet) {
1401       return 2; // all detectors in road are bad
1402     } else {
1403       return 3; // at least one is bad
1404     }
1405   }
1406
1407   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
1408   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
1409   if (!nChipsInRoad) return 0;
1410
1411   Bool_t anyBad=kFALSE,anyGood=kFALSE;
1412   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
1413     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
1414     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
1415     if (det.IsChipBad(chipsInRoad[iCh])) {
1416       anyBad=kTRUE;
1417     } else {
1418       anyGood=kTRUE;
1419     }
1420   }
1421
1422   if (!anyGood) {
1423     if(!touchNeighbourDet) {
1424       AliDebug(2,"all bad in road");
1425       return 2;  // all chips in road are bad
1426     } else {
1427       return 3; // at least a bad chip in road
1428     }
1429   }
1430
1431   if (anyBad) {
1432     AliDebug(2,"at least a bad in road");
1433     return 3; // at least a bad chip in road
1434   }
1435
1436
1437   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
1438       || ilayer==4 || ilayer==5     // SSD
1439       || !noClusters) return 0;
1440
1441   // There are no clusters in road: check if there is at least
1442   // a bad SPD pixel or SDD anode
1443
1444   Int_t idetInITS=idet;
1445   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
1446
1447   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
1448     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
1449     return 3;
1450   }
1451   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
1452
1453   return 0;
1454 }
1455 //------------------------------------------------------------------------
1456 Bool_t AliITStrackerHLT::LocalModuleCoord(Int_t ilayer,Int_t idet,
1457                                        const AliHLTITSTrack *track,
1458                                        Float_t &xloc,Float_t &zloc) const {
1459   //-----------------------------------------------------------------
1460   // Gives position of track in local module ref. frame
1461   //-----------------------------------------------------------------
1462
1463   xloc=0.;
1464   zloc=0.;
1465
1466   if(idet<0) return kFALSE;
1467
1468   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
1469
1470   Int_t lad = Int_t(idet/ndet) + 1;
1471
1472   Int_t det = idet - (lad-1)*ndet + 1;
1473
1474   Double_t xyzGlob[3],xyzLoc[3];
1475
1476   AliHLTITSDetector &detector = fgLayers[ilayer].GetDetector(idet);
1477   // take into account the misalignment: xyz at real detector plane
1478   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
1479
1480   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
1481
1482   xloc = (Float_t)xyzLoc[0];
1483   zloc = (Float_t)xyzLoc[2];
1484
1485   return kTRUE;
1486 }
1487
1488 //------------------------------------------------------------------------
1489 Bool_t AliITStrackerHLT::ComputeRoad(AliHLTITSTrack* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1490   //--------------------------------------------------------------------
1491   // This function computes the rectangular road for this track
1492   //--------------------------------------------------------------------
1493
1494
1495   AliHLTITSDetector &det = fgLayers[ilayer].GetDetector(idet);
1496   // take into account the misalignment: propagate track to misaligned detector plane
1497   if (!TransportToPhiX( track, det.GetPhi(),det.GetRmisal() ) ) return kFALSE;
1498
1499   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1500                     TMath::Sqrt(track->GetSigmaZ2() +
1501                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1502                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1503                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1504   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1505                     TMath::Sqrt(track->GetSigmaY2() +
1506                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1507                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1508                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1509
1510   // track at boundary between detectors, enlarge road
1511   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1512   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1513        (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1514        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1515        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1516     Float_t tgl = TMath::Abs(track->GetTgl());
1517     if (tgl > 1.) tgl=1.;
1518     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1519     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1520     Float_t snp = TMath::Abs(track->GetSnp());
1521     if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1522     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1523   } // boundary
1524
1525   // add to the road a term (up to 2-3 mm) to deal with misalignments
1526   dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1527   dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1528
1529   Double_t r = fgLayers[ilayer].GetR();
1530   zmin = track->GetZ() - dz;
1531   zmax = track->GetZ() + dz;
1532   ymin = track->GetY() + r*det.GetPhi() - dy;
1533   ymax = track->GetY() + r*det.GetPhi() + dy;
1534
1535   // bring track back to idead detector plane
1536   if (!TransportToPhiX( track, det.GetPhi(),det.GetR())) return kFALSE;
1537
1538   return kTRUE;
1539 }