ITS/AliITStrackerMI.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$ */
  16
  17 //-------------------------------------------------------------------------
  18 //               Implementation of the ITS tracker class
  19 //    It reads AliITSRecPoint clusters and creates AliITStrackMI 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 #include <TVector3.h>
  36 #include <TBits.h>
  37
  38 #include "AliLog.h"
  39 #include "AliGeomManager.h"
  40 #include "AliITSPlaneEff.h"
  41 #include "AliTrackPointArray.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDV0Params.h"
  44 #include "AliESDtrack.h"
  45 #include "AliV0.h"
  46 #include "AliITSChannelStatus.h"
  47 #include "AliITSDetTypeRec.h"
  48 #include "AliITSRecPoint.h"
  49 #include "AliITSRecPointContainer.h"
  50 #include "AliITSgeomTGeo.h"
  51 #include "AliITSReconstructor.h"
  52 #include "AliITSClusterParam.h"
  53 #include "AliITSsegmentation.h"
  54 #include "AliITSCalibration.h"
  55 #include "AliITSPlaneEffSPD.h"
  56 #include "AliITSPlaneEffSDD.h"
  57 #include "AliITSPlaneEffSSD.h"
  58 #include "AliITSV0Finder.h"
  59 #include "AliITStrackerMI.h"
  60 #include "AliMathBase.h"
  61
  62
  63 ClassImp(AliITStrackerMI)
  64
  65 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
  66
  67 AliITStrackerMI::AliITStrackerMI():AliTracker(),
  68 fI(0),
  69 fBestTrack(),
  70 fTrackToFollow(),
  71 fTrackHypothesys(),
  72 fBestHypothesys(),
  73 fOriginal(),
  74 fCurrentEsdTrack(),
  75 fPass(0),
  76 fAfterV0(kFALSE),
  77 fLastLayerToTrackTo(0),
  78 fCoefficients(0),
  79 fEsd(0),
  80 fTrackingPhase("Default"),
  81 fUseTGeo(3),
  82 fNtracks(0),
  83 fFlagFakes(kFALSE),
  84 fSelectBestMIP03(kFALSE),
  85 fUseImproveKalman(kFALSE),
  86 fxOverX0Pipe(-1.),
  87 fxTimesRhoPipe(-1.),
  88 fxOverX0PipeTrks(0),
  89 fxTimesRhoPipeTrks(0),
  90 fxOverX0ShieldTrks(0),
  91 fxTimesRhoShieldTrks(0),
  92 fxOverX0LayerTrks(0),
  93 fxTimesRhoLayerTrks(0),
  94 fDebugStreamer(0),
  95 fITSChannelStatus(0),
  96 fkDetTypeRec(0),
  97 fPlaneEff(0),
  98 fSPDChipIntPlaneEff(0),
  99 fITSPid(0)
 100
 101  {
 102   //Default constructor
 103   Int_t i;
 104   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
 105   for(i=0;i<2;i++) {
 106     fxOverX0Shield[i]=-1.;
 107     fxTimesRhoShield[i]=-1.;
 108     fConstraint[i]=0;
 109   }
 110   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 111   fOriginal.SetOwner();
 112   for(i=0;i<AliITSgeomTGeo::kNLayers;i++)fForceSkippingOfLayer[i]=0;
 113   for(i=0;i<100000;i++)fBestTrackIndex[i]=0;
 114   fITSPid=new AliITSPIDResponse();
 115
 116 }
 117 //------------------------------------------------------------------------
 118 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
 119 fI(AliITSgeomTGeo::GetNLayers()),
 120 fBestTrack(),
 121 fTrackToFollow(),
 122 fTrackHypothesys(),
 123 fBestHypothesys(),
 124 fOriginal(),
 125 fCurrentEsdTrack(),
 126 fPass(0),
 127 fAfterV0(kFALSE),
 128 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
 129 fCoefficients(0),
 130 fEsd(0),
 131 fTrackingPhase("Default"),
 132 fUseTGeo(3),
 133 fNtracks(0),
 134 fFlagFakes(kFALSE),
 135 fSelectBestMIP03(kFALSE),
 136 fUseImproveKalman(kFALSE),
 137 fxOverX0Pipe(-1.),
 138 fxTimesRhoPipe(-1.),
 139 fxOverX0PipeTrks(0),
 140 fxTimesRhoPipeTrks(0),
 141 fxOverX0ShieldTrks(0),
 142 fxTimesRhoShieldTrks(0),
 143 fxOverX0LayerTrks(0),
 144 fxTimesRhoLayerTrks(0),
 145 fDebugStreamer(0),
 146 fITSChannelStatus(0),
 147 fkDetTypeRec(0),
 148 fPlaneEff(0),
 149 fSPDChipIntPlaneEff(0),
 150 fITSPid(0) {
 151   //--------------------------------------------------------------------
 152   //This is the AliITStrackerMI constructor
 153   //--------------------------------------------------------------------
 154   if (geom) {
 155     AliWarning("\"geom\" is actually a dummy argument !");
 156   }
 157
 158   fOriginal.SetOwner();
 159   fCoefficients = 0;
 160   fAfterV0     = kFALSE;
 161
 162   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 163     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
 164     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
 165
 166     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
 167     AliITSgeomTGeo::GetTranslation(i,1,1,xyz);
 168     Double_t poff=TMath::ATan2(y,x);
 169     Double_t zoff=z;
 170
 171     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 172     Double_t r=TMath::Sqrt(x*x + y*y);
 173
 174     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 175     r += TMath::Sqrt(x*x + y*y);
 176     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 177     r += TMath::Sqrt(x*x + y*y);
 178     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 179     r += TMath::Sqrt(x*x + y*y);
 180     r*=0.25;
 181
 182     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
 183
 184     for (Int_t j=1; j<nlad+1; j++) {
 185       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 186         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 187         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 188         m.Multiply(tm);
 189         Double_t txyz[3]={0.};
 190         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 191         m.LocalToMaster(txyz,xyz);
 192         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 193         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 194
 195         if (phi<0) phi+=TMath::TwoPi();
 196         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 197
 198         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 199         new(&det) AliITSdetector(r,phi);
 200         // compute the real radius (with misalignment)
 201         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
 202         mmisal.Multiply(tm);
 203         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 204         mmisal.LocalToMaster(txyz,xyz);
 205         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 206         det.SetRmisal(rmisal);
 207
 208       } // end loop on detectors
 209     } // end loop on ladders
 210     fForceSkippingOfLayer[i-1] = 0;
 211   } // end loop on layers
 212
 213
 214   fI=AliITSgeomTGeo::GetNLayers();
 215
 216   fPass=0;
 217   fConstraint[0]=1; fConstraint[1]=0;
 218
 219   Double_t xyzVtx[]={AliITSReconstructor::GetRecoParam()->GetXVdef(),
 220                      AliITSReconstructor::GetRecoParam()->GetYVdef(),
 221                      AliITSReconstructor::GetRecoParam()->GetZVdef()};
 222   Double_t ersVtx[]={AliITSReconstructor::GetRecoParam()->GetSigmaXVdef(),
 223                      AliITSReconstructor::GetRecoParam()->GetSigmaYVdef(),
 224                      AliITSReconstructor::GetRecoParam()->GetSigmaZVdef()};
 225   SetVertex(xyzVtx,ersVtx);
 226
 227   fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
 228   for (Int_t i=0;i<100000;i++){
 229     fBestTrackIndex[i]=0;
 230   }
 231
 232   // store positions of centre of SPD modules (in z)
 233   //  The convetion is that fSPDdetzcentre is ordered from -z to +z
 234   Double_t tr[3];
 235   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 236   if (tr[2]<0) { // old geom (up to v5asymmPPR)
 237     AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 238     fSPDdetzcentre[0] = tr[2];
 239     AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 240     fSPDdetzcentre[1] = tr[2];
 241     AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 242     fSPDdetzcentre[2] = tr[2];
 243     AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 244     fSPDdetzcentre[3] = tr[2];
 245   } else { // new geom (from v11Hybrid)
 246     AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 247     fSPDdetzcentre[0] = tr[2];
 248     AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 249     fSPDdetzcentre[1] = tr[2];
 250     AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 251     fSPDdetzcentre[2] = tr[2];
 252     AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 253     fSPDdetzcentre[3] = tr[2];
 254   }
 255
 256   fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
 257   if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
 258     AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
 259     fUseTGeo = 3;
 260   }
 261
 262   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 263   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 264
 265   if (AliITSReconstructor::GetRecoParam()->GetESDV0Params()->StreamLevel()>0)
 266     fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 267
 268   // only for plane efficiency evaluation
 269   if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 270       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
 271     Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
 272     if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
 273       AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
 274     if (iplane<2) {
 275       fPlaneEff = new AliITSPlaneEffSPD();
 276       fSPDChipIntPlaneEff = new Bool_t[AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip];
 277       for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
 278     }
 279     else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
 280     else fPlaneEff = new AliITSPlaneEffSSD();
 281     if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
 282        if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
 283     if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
 284   }
 285   //
 286   // RS
 287   fSelectBestMIP03  = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
 288   fFlagFakes        = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
 289   fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
 290   //
 291   fITSPid=new AliITSPIDResponse();
 292 }
 293 /*
 294 //------------------------------------------------------------------------
 295 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
 296 fI(tracker.fI),
 297 fBestTrack(tracker.fBestTrack),
 298 fTrackToFollow(tracker.fTrackToFollow),
 299 fTrackHypothesys(tracker.fTrackHypothesys),
 300 fBestHypothesys(tracker.fBestHypothesys),
 301 fOriginal(tracker.fOriginal),
 302 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
 303 fPass(tracker.fPass),
 304 fAfterV0(tracker.fAfterV0),
 305 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
 306 fCoefficients(tracker.fCoefficients),
 307 fEsd(tracker.fEsd),
 308 fTrackingPhase(tracker.fTrackingPhase),
 309 fUseTGeo(tracker.fUseTGeo),
 310 fNtracks(tracker.fNtracks),
 311 fFlagFakes(tracker.fFlagFakes),
 312 fSelectBestMIP03(tracker.fSelectBestMIP03),
 313 fxOverX0Pipe(tracker.fxOverX0Pipe),
 314 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
 315 fxOverX0PipeTrks(0),
 316 fxTimesRhoPipeTrks(0),
 317 fxOverX0ShieldTrks(0),
 318 fxTimesRhoShieldTrks(0),
 319 fxOverX0LayerTrks(0),
 320 fxTimesRhoLayerTrks(0),
 321 fDebugStreamer(tracker.fDebugStreamer),
 322 fITSChannelStatus(tracker.fITSChannelStatus),
 323 fkDetTypeRec(tracker.fkDetTypeRec),
 324 fPlaneEff(tracker.fPlaneEff) {
 325   //Copy constructor
 326   fOriginal.SetOwner();
 327   Int_t i;
 328   for(i=0;i<4;i++) {
 329     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
 330   }
 331   for(i=0;i<6;i++) {
 332     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
 333     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
 334   }
 335   for(i=0;i<2;i++) {
 336     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
 337     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
 338   }
 339 }
 340
 341 //------------------------------------------------------------------------
 342 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
 343   //Assignment operator
 344   this->~AliITStrackerMI();
 345   new(this) AliITStrackerMI(tracker);
 346   return *this;
 347 }
 348 */
 349 //------------------------------------------------------------------------
 350 AliITStrackerMI::~AliITStrackerMI()
 351 {
 352   //
 353   //destructor
 354   //
 355   if (fCoefficients) delete [] fCoefficients;
 356   DeleteTrksMaterialLUT();
 357   if (fDebugStreamer) {
 358     //fDebugStreamer->Close();
 359     delete fDebugStreamer;
 360   }
 361   if(fITSChannelStatus) delete fITSChannelStatus;
 362   if(fPlaneEff) delete fPlaneEff;
 363   if(fITSPid) delete fITSPid;
 364   if (fSPDChipIntPlaneEff) delete [] fSPDChipIntPlaneEff;
 365
 366 }
 367 //------------------------------------------------------------------------
 368 void AliITStrackerMI::ReadBadFromDetTypeRec() {
 369   //--------------------------------------------------------------------
 370   //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
 371   //i.e. from OCDB
 372   //--------------------------------------------------------------------
 373
 374   if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
 375
 376   Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
 377
 378   if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
 379
 380   // ITS channels map
 381   if(fITSChannelStatus) delete fITSChannelStatus;
 382   fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
 383
 384   // ITS detectors and chips
 385   Int_t i=0,j=0,k=0,ndet=0;
 386   for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 387     Int_t nBadDetsPerLayer=0;
 388     ndet=AliITSgeomTGeo::GetNDetectors(i);
 389     for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
 390       for (k=1; k<ndet+1; k++) {
 391         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 392         det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
 393         if(det.IsBad()) {nBadDetsPerLayer++;}
 394       } // end loop on detectors
 395     } // end loop on ladders
 396     AliInfo(Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
 397   } // end loop on layers
 398
 399   return;
 400 }
 401 //------------------------------------------------------------------------
 402 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
 403   //--------------------------------------------------------------------
 404   //This function loads ITS clusters
 405   //--------------------------------------------------------------------
 406
 407   TClonesArray *clusters = NULL;
 408   AliITSRecPointContainer* rpcont=AliITSRecPointContainer::Instance();
 409   clusters=rpcont->FetchClusters(0,cTree);
 410   if(!clusters) return 1;
 411
 412   if(!(rpcont->IsSPDActive() || rpcont->IsSDDActive() || rpcont->IsSSDActive())){
 413       AliError("ITS is not in a known running configuration: SPD, SDD and SSD are not active");
 414       return 1;
 415   }
 416   Int_t i=0,j=0,ndet=0;
 417   Int_t detector=0;
 418   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 419     ndet=fgLayers[i].GetNdetectors();
 420     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 421     for (; j<jmax; j++) {
 422       //      if (!cTree->GetEvent(j)) continue;
 423       clusters = rpcont->UncheckedGetClusters(j);
 424       if(!clusters)continue;
 425       Int_t ncl=clusters->GetEntriesFast();
 426       SignDeltas(clusters,GetZ());
 427
 428       while (ncl--) {
 429         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 430         detector=c->GetDetectorIndex();
 431
 432         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 433
 434         Int_t retval = fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 435         if(retval) {
 436           AliWarning(Form("Too many clusters on layer %d!",i));
 437           break;
 438         }
 439       }
 440
 441       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 442       // zwindow cm from the dead zone
 443       //  This method assumes that fSPDdetzcentre is ordered from -z to +z
 444       if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
 445         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
 446           Int_t lab[4]   = {0,0,0,detector};
 447           Int_t info[3]  = {0,0,i};
 448           Float_t q      = 0.; // this identifies virtual clusters
 449           Float_t hit[6] = {xdead,
 450                             0.,
 451                             AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
 452                             AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
 453                             q,
 454                             0.};
 455           Bool_t local   = kTRUE;
 456           Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
 457           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 458           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 459             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 460           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 461           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 462             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 463           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 464           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 465             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 466           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 467           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 468             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 469           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 470           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 471             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 472           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
 473           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 474             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 475         }
 476       } // "virtual" clusters in SPD
 477
 478     }
 479     //
 480     fgLayers[i].ResetRoad(); //road defined by the cluster density
 481     fgLayers[i].SortClusters();
 482   }
 483
 484   // check whether we have to skip some layers
 485   SetForceSkippingOfLayer();
 486
 487   return 0;
 488 }
 489 //------------------------------------------------------------------------
 490 void AliITStrackerMI::UnloadClusters() {
 491   //--------------------------------------------------------------------
 492   //This function unloads ITS clusters
 493   //--------------------------------------------------------------------
 494   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
 495 }
 496 //------------------------------------------------------------------------
 497 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
 498   //--------------------------------------------------------------------
 499   // Publishes all pointers to clusters known to the tracker into the
 500   // passed object array.
 501   // The ownership is not transfered - the caller is not expected to delete
 502   // the clusters.
 503   //--------------------------------------------------------------------
 504
 505   for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 506     for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
 507       AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
 508       array->AddLast(cl);
 509     }
 510   }
 511
 512   return;
 513 }
 514 //------------------------------------------------------------------------
 515 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
 516   //--------------------------------------------------------------------
 517   // Correction for the material between the TPC and the ITS
 518   //--------------------------------------------------------------------
 519   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction
 520       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
 521       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
 522       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
 523   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
 524       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
 525       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
 526       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
 527   } else {
 528     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
 529     return 0;
 530   }
 531
 532   return 1;
 533 }
 534 //------------------------------------------------------------------------
 535 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
 536   //--------------------------------------------------------------------
 537   // This functions reconstructs ITS tracks
 538   // The clusters must be already loaded !
 539   //--------------------------------------------------------------------
 540
 541   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
 542
 543   fTrackingPhase="Clusters2Tracks";
 544   //
 545   // RS
 546   fSelectBestMIP03  = kFALSE;//AliITSReconstructor::GetRecoParam()->GetSelectBestMIP03();
 547   fFlagFakes        = AliITSReconstructor::GetRecoParam()->GetFlagFakes();
 548   fUseImproveKalman = AliITSReconstructor::GetRecoParam()->GetUseImproveKalman();
 549   //
 550   TObjArray itsTracks(15000);
 551   fOriginal.Clear();
 552   fEsd = event;         // store pointer to the esd
 553
 554   // temporary (for cosmics)
 555   if(event->GetVertex()) {
 556     TString title = event->GetVertex()->GetTitle();
 557     if(title.Contains("cosmics")) {
 558       Double_t xyz[3]={GetX(),GetY(),GetZ()};
 559       Double_t exyz[3]={0.1,0.1,0.1};
 560       SetVertex(xyz,exyz);
 561     }
 562   }
 563   // temporary
 564   Int_t noesd = 0;
 565   {/* Read ESD tracks */
 566     Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 567     Int_t nentr=event->GetNumberOfTracks();
 568     noesd=nentr;
 569     //    Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 570     while (nentr--) {
 571       AliESDtrack *esd=event->GetTrack(nentr);
 572       //  ---- for debugging:
 573       //if(TMath::Abs(esd->GetX()-83.65)<0.1) { FILE *f=fopen("tpc.dat","a"); fprintf(f,"%f %f %f %f %f %f\n",(Float_t)event->GetEventNumberInFile(),(Float_t)TMath::Abs(esd->GetLabel()),(Float_t)esd->GetX(),(Float_t)esd->GetY(),(Float_t)esd->GetZ(),(Float_t)esd->Pt()); fclose(f); }
 574
 575       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 576       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 577       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 578       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
 579       AliITStrackMI *t = new AliITStrackMI(*esd);
 580       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
 581       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
 582
 583
 584       // look at the ESD mass hypothesys !
 585       if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
 586       // write expected q
 587       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 588
 589       if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
 590         //track - can be  V0 according to TPC
 591       } else {
 592         if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
 593           delete t;
 594           continue;
 595         }
 596         if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
 597           delete t;
 598           continue;
 599         }
 600         if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
 601           delete t;
 602           continue;
 603         }
 604         if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 605           delete t;
 606           continue;
 607         }
 608       }
 609       t->SetReconstructed(kFALSE);
 610       itsTracks.AddLast(t);
 611       fOriginal.AddLast(t);
 612     }
 613   } /* End Read ESD tracks */
 614
 615   itsTracks.Sort();
 616   fOriginal.Sort();
 617   Int_t nentr=itsTracks.GetEntriesFast();
 618   fTrackHypothesys.Expand(nentr);
 619   fBestHypothesys.Expand(nentr);
 620   MakeCoefficients(nentr);
 621   if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
 622   Int_t ntrk=0;
 623   // THE TWO TRACKING PASSES
 624   for (fPass=0; fPass<2; fPass++) {
 625      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 626      for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
 627        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
 628        if (t==0) continue;              //this track has been already tracked
 629        //cout<<"========== "<<fPass<<"    "<<fCurrentEsdTrack<<" =========\n";
 630        if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue;  //this track was  already  "succesfully" reconstructed
 631        Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz);              //I.B.
 632        if (fConstraint[fPass]) {
 633          if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
 634              TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
 635        }
 636
 637        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 638        AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
 639        fI = 6;
 640        ResetTrackToFollow(*t);
 641        ResetBestTrack();
 642
 643        FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
 644
 645
 646        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
 647        //
 648        AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
 649        if (!besttrack) continue;
 650        besttrack->SetLabel(tpcLabel);
 651        //       besttrack->CookdEdx();
 652        CookdEdx(besttrack);
 653        besttrack->SetFakeRatio(1.);
 654        CookLabel(besttrack,0.); //For comparison only
 655        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
 656        t->SetWinner(besttrack);
 657
 658        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain
 659
 660        t->SetReconstructed(kTRUE);
 661        ntrk++;
 662        AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
 663      }
 664      GetBestHypothesysMIP(itsTracks);
 665   } // end loop on the two tracking passes
 666   //
 667   if (fFlagFakes) FlagFakes(itsTracks);
 668   //
 669   if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
 670   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
 671   fAfterV0 = kTRUE;
 672   //
 673   itsTracks.Clear();
 674   //
 675   Int_t entries = fTrackHypothesys.GetEntriesFast();
 676   for (Int_t ientry=0; ientry<entries; ientry++) {
 677     TObjArray * array =(TObjArray*)fTrackHypothesys.At(ientry);
 678     if (array) array->Delete();
 679     delete fTrackHypothesys.RemoveAt(ientry);
 680   }
 681
 682   fTrackHypothesys.Delete();
 683   entries = fBestHypothesys.GetEntriesFast();
 684   for (Int_t ientry=0; ientry<entries; ientry++) {
 685     TObjArray * array =(TObjArray*)fBestHypothesys.At(ientry);
 686     if (array) array->Delete();
 687     delete fBestHypothesys.RemoveAt(ientry);
 688   }
 689   fBestHypothesys.Delete();
 690   fOriginal.Clear();
 691   delete [] fCoefficients;
 692   fCoefficients=0;
 693   DeleteTrksMaterialLUT();
 694
 695   AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
 696
 697   fTrackingPhase="Default";
 698
 699   return 0;
 700 }
 701 //------------------------------------------------------------------------
 702 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
 703   //--------------------------------------------------------------------
 704   // This functions propagates reconstructed ITS tracks back
 705   // The clusters must be loaded !
 706   //--------------------------------------------------------------------
 707   fTrackingPhase="PropagateBack";
 708   Int_t nentr=event->GetNumberOfTracks();
 709   //  Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 710   double bz0 = GetBz();
 711   const double kWatchStep=10.; // for larger steps watch arc vs segment difference
 712   //
 713   Int_t ntrk=0;
 714   for (Int_t i=0; i<nentr; i++) {
 715      AliESDtrack *esd=event->GetTrack(i);
 716
 717      // Start time integral and add distance from current position to vertex
 718      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 719      AliITStrackMI t(*esd);
 720      Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
 721      t.GetXYZ(xyzTrk);
 722      Double_t dst2 = 0.;
 723      {
 724        double dxs = xyzTrk[0] - xyzVtx[0];
 725        double dys = xyzTrk[1] - xyzVtx[1];
 726        double dzs = xyzTrk[2] - xyzVtx[2];
 727        // RS: for large segment steps d use approximation of cicrular arc s by
 728        // s = 2R*asin(d/2R) = d/p asin(p) \approx d/p (p + 1/6 p^3) = d (1+1/6 p^2)
 729        // where R is the track radius, p = d/2R = 2C*d (C is the abs curvature)
 730        // Hence s^2/d^2 = (1+1/6 p^2)^2
 731        dst2 = dxs*dxs + dys*dys;
 732        if (dst2 > kWatchStep*kWatchStep) { // correct circular part for arc/segment factor
 733          double crv = TMath::Abs(esd->GetC(bz0));
 734          double fcarc = 1.+crv*crv*dst2/6.;
 735          dst2 *= fcarc*fcarc;
 736        }
 737        dst2 += dzs*dzs;
 738      }
 739      t.StartTimeIntegral();
 740      t.AddTimeStep(TMath::Sqrt(dst2));
 741      //
 742      // transfer the time integral to ESD track
 743      esd->SetStatus(AliESDtrack::kTIME);
 744      Double_t times[10];t.GetIntegratedTimes(times); esd->SetIntegratedTimes(times);
 745      esd->SetIntegratedLength(t.GetIntegratedLength());
 746      //
 747      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 748
 749      t.SetExpQ(TMath::Max(0.8*t.GetESDtrack()->GetTPCsignal(),30.));
 750      ResetTrackToFollow(t);
 751      //
 752      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
 753      if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,&t)) {
 754        if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) continue;
 755        fTrackToFollow.SetLabel(t.GetLabel());
 756        //fTrackToFollow.CookdEdx();
 757        CookLabel(&fTrackToFollow,0.); //For comparison only
 758        fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 759        //UseClusters(&fTrackToFollow);
 760        ntrk++;
 761      }
 762   }
 763
 764   AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
 765
 766   fTrackingPhase="Default";
 767
 768   return 0;
 769 }
 770 //------------------------------------------------------------------------
 771 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
 772   //--------------------------------------------------------------------
 773   // This functions refits ITS tracks using the
 774   // "inward propagated" TPC tracks
 775   // The clusters must be loaded !
 776   //--------------------------------------------------------------------
 777   fTrackingPhase="RefitInward";
 778
 779   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
 780
 781   Bool_t doExtra=AliITSReconstructor::GetRecoParam()->GetSearchForExtraClusters();
 782   if(!doExtra) AliDebug(2,"Do not search for extra clusters");
 783
 784   Int_t nentr=event->GetNumberOfTracks();
 785   //  Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 786
 787   // only for PlaneEff and in case of SPD (for FO studies)
 788   if( AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 789       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0 &&
 790       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()<2) {
 791       for (UInt_t i=0; i<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip; i++) fSPDChipIntPlaneEff[i]=kFALSE;
 792   }
 793
 794   Int_t ntrk=0;
 795   for (Int_t i=0; i<nentr; i++) {
 796     AliESDtrack *esd=event->GetTrack(i);
 797
 798     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 799     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 800     if (esd->GetStatus()&AliESDtrack::kTPCout)
 801       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 802
 803     AliITStrackMI *t = new AliITStrackMI(*esd);
 804
 805     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 806     if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 807        delete t;
 808        continue;
 809     }
 810
 811     ResetTrackToFollow(*t);
 812     fTrackToFollow.ResetClusters();
 813
 814     // ITS standalone tracks
 815     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) {
 816       fTrackToFollow.ResetCovariance(10.);
 817       // protection for loopers that can have parameters screwed up
 818       if(TMath::Abs(fTrackToFollow.GetY())>1000. ||
 819          TMath::Abs(fTrackToFollow.GetZ())>1000.) {
 820         delete t;
 821         continue;
 822       }
 823     }
 824
 825     //Refitting...
 826     Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 827                AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
 828
 829     AliDebug(2,Form("Refit LABEL %d  %d",t->GetLabel(),t->GetNumberOfClusters()));
 830     if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,doExtra,pe)) {
 831        AliDebug(2,"  refit OK");
 832        fTrackToFollow.SetLabel(t->GetLabel());
 833        //       fTrackToFollow.CookdEdx();
 834        CookdEdx(&fTrackToFollow);
 835
 836        CookLabel(&fTrackToFollow,0.0); //For comparison only
 837
 838        //The beam pipe
 839        if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 840          fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
 841          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
 842          //printf("                                       %d\n",esdTrack->GetITSModuleIndex(0));
 843          //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
 844          Double_t r[3]={0.,0.,0.};
 845          Double_t maxD=3.;
 846          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
 847          ntrk++;
 848        }
 849     }
 850     delete t;
 851   }
 852
 853   AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
 854
 855   fTrackingPhase="Default";
 856
 857   return 0;
 858 }
 859 //------------------------------------------------------------------------
 860 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 861   //--------------------------------------------------------------------
 862   //       Return pointer to a given cluster
 863   //--------------------------------------------------------------------
 864   Int_t l=(index & 0xf0000000) >> 28;
 865   Int_t c=(index & 0x0fffffff) >> 00;
 866   return fgLayers[l].GetCluster(c);
 867 }
 868 //------------------------------------------------------------------------
 869 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 870   //--------------------------------------------------------------------
 871   // Get track space point with index i
 872   //--------------------------------------------------------------------
 873
 874   Int_t l=(index & 0xf0000000) >> 28;
 875   Int_t c=(index & 0x0fffffff) >> 00;
 876   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 877   Int_t idet = cl->GetDetectorIndex();
 878
 879   Float_t xyz[3];
 880   Float_t cov[6];
 881   cl->GetGlobalXYZ(xyz);
 882   cl->GetGlobalCov(cov);
 883   p.SetXYZ(xyz, cov);
 884   p.SetCharge(cl->GetQ());
 885   p.SetDriftTime(cl->GetDriftTime());
 886   p.SetChargeRatio(cl->GetChargeRatio());
 887   p.SetClusterType(cl->GetClusterType());
 888   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 889   switch (l) {
 890   case 0:
 891     iLayer = AliGeomManager::kSPD1;
 892     break;
 893   case 1:
 894     iLayer = AliGeomManager::kSPD2;
 895     break;
 896   case 2:
 897     iLayer = AliGeomManager::kSDD1;
 898     break;
 899   case 3:
 900     iLayer = AliGeomManager::kSDD2;
 901     break;
 902   case 4:
 903     iLayer = AliGeomManager::kSSD1;
 904     break;
 905   case 5:
 906     iLayer = AliGeomManager::kSSD2;
 907     break;
 908   default:
 909     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 910     break;
 911   };
 912   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 913   p.SetVolumeID((UShort_t)volid);
 914   return kTRUE;
 915 }
 916 //------------------------------------------------------------------------
 917 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
 918                         AliTrackPoint& p, const AliESDtrack *t) {
 919   //--------------------------------------------------------------------
 920   // Get track space point with index i
 921   // (assign error estimated during the tracking)
 922   //--------------------------------------------------------------------
 923
 924   Int_t l=(index & 0xf0000000) >> 28;
 925   Int_t c=(index & 0x0fffffff) >> 00;
 926   const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 927   Int_t idet = cl->GetDetectorIndex();
 928
 929   const AliITSdetector &det=fgLayers[l].GetDetector(idet);
 930
 931   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
 932   Float_t detxy[2];
 933   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
 934   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
 935   Double_t alpha = t->GetAlpha();
 936   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
 937   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe+cl->GetX(),GetBz()));
 938   phi += alpha-det.GetPhi();
 939   Float_t tgphi = TMath::Tan(phi);
 940
 941   Float_t tgl = t->GetTgl(); // tgl about const along track
 942   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
 943
 944   Float_t errtrky,errtrkz,covyz;
 945   Bool_t addMisalErr=kFALSE;
 946   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
 947
 948   Float_t xyz[3];
 949   Float_t cov[6];
 950   cl->GetGlobalXYZ(xyz);
 951   //  cl->GetGlobalCov(cov);
 952   Float_t pos[3] = {0.,0.,0.};
 953   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
 954   tmpcl.GetGlobalCov(cov);
 955
 956   p.SetXYZ(xyz, cov);
 957   p.SetCharge(cl->GetQ());
 958   p.SetDriftTime(cl->GetDriftTime());
 959   p.SetChargeRatio(cl->GetChargeRatio());
 960   p.SetClusterType(cl->GetClusterType());
 961
 962   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 963   switch (l) {
 964   case 0:
 965     iLayer = AliGeomManager::kSPD1;
 966     break;
 967   case 1:
 968     iLayer = AliGeomManager::kSPD2;
 969     break;
 970   case 2:
 971     iLayer = AliGeomManager::kSDD1;
 972     break;
 973   case 3:
 974     iLayer = AliGeomManager::kSDD2;
 975     break;
 976   case 4:
 977     iLayer = AliGeomManager::kSSD1;
 978     break;
 979   case 5:
 980     iLayer = AliGeomManager::kSSD2;
 981     break;
 982   default:
 983     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 984     break;
 985   };
 986   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 987
 988   p.SetVolumeID((UShort_t)volid);
 989   return kTRUE;
 990 }
 991 //------------------------------------------------------------------------
 992 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
 993 {
 994   //--------------------------------------------------------------------
 995   // Follow prolongation tree
 996   //--------------------------------------------------------------------
 997   //
 998   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
 999   Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
1000
1001
1002   AliESDtrack * esd = otrack->GetESDtrack();
1003   if (esd->GetV0Index(0)>0) {
1004     // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
1005     //                      mapping of ESD track is different as ITS track in Containers
1006     //                      Need something more stable
1007     //                      Indexes are set back again to the ESD track indexes in UpdateTPCV0
1008     for (Int_t i=0;i<3;i++){
1009       Int_t  index = esd->GetV0Index(i);
1010       if (index==0) break;
1011       AliESDv0 * vertex = fEsd->GetV0(index);
1012       if (vertex->GetStatus()<0) continue;     // rejected V0
1013       //
1014       if (esd->GetSign()>0) {
1015         vertex->SetIndex(0,esdindex);
1016       } else {
1017         vertex->SetIndex(1,esdindex);
1018       }
1019     }
1020   }
1021   TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
1022   if (!bestarray){
1023     bestarray = new TObjArray(5);
1024     bestarray->SetOwner();
1025     fBestHypothesys.AddAt(bestarray,esdindex);
1026   }
1027
1028   //
1029   //setup tree of the prolongations
1030   //
1031   const int kMaxTr = 100; //RS
1032   static AliITStrackMI tracks[7][kMaxTr];
1033   AliITStrackMI *currenttrack;
1034   static AliITStrackMI currenttrack1;
1035   static AliITStrackMI currenttrack2;
1036   static AliITStrackMI backuptrack;
1037   Int_t ntracks[7];
1038   Int_t nindexes[7][kMaxTr];
1039   Float_t normalizedchi2[kMaxTr];
1040   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
1041   otrack->SetNSkipped(0);
1042   new (&(tracks[6][0])) AliITStrackMI(*otrack);
1043   ntracks[6]=1;
1044   for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
1045   Int_t modstatus = 1; // found
1046   Float_t xloc,zloc;
1047   //
1048   //
1049   // follow prolongations
1050   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
1051     AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
1052     fI = ilayer;
1053     //
1054     AliITSlayer &layer=fgLayers[ilayer];
1055     Double_t r = layer.GetR();
1056     ntracks[ilayer]=0;
1057     //
1058     //
1059     Int_t nskipped=0;
1060     Float_t nused =0;
1061     for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
1062       //      printf("LR %d Tr:%d NSeeds: %d\n",ilayer,itrack,ntracks[ilayer+1]);
1063       //set current track
1064       if (ntracks[ilayer]>=kMaxTr) break;
1065       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
1066       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
1067       if (ntracks[ilayer]>15+ilayer){
1068         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
1069         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
1070       }
1071
1072       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
1073
1074       // material between SSD and SDD, SDD and SPD
1075       if (ilayer==3)
1076         if(!CorrectForShieldMaterial(&currenttrack1,"SDD","inward")) continue;
1077       if (ilayer==1)
1078         if(!CorrectForShieldMaterial(&currenttrack1,"SPD","inward")) continue;
1079
1080       // detector number
1081       Double_t phi,z;
1082       if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
1083       Int_t idet=layer.FindDetectorIndex(phi,z);
1084
1085       Double_t trackGlobXYZ1[3];
1086       if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1087
1088       // Get the budget to the primary vertex for the current track being prolonged
1089       Double_t budgetToPrimVertex = 0;
1090       double xMSLrs[9],x2X0MSLrs[9]; // needed for ImproveKalman
1091       int nMSLrs = 0;
1092       //
1093       if (fUseImproveKalman) nMSLrs = GetEffectiveThicknessLbyL(xMSLrs,x2X0MSLrs);
1094       else budgetToPrimVertex = GetEffectiveThickness();
1095       //
1096       // check if we allow a prolongation without point
1097       Int_t skip = CheckSkipLayer(&currenttrack1,ilayer,idet);
1098       if (skip) {
1099         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1100         // propagate to the layer radius
1101         Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1102         if(!vtrack->Propagate(xToGo)) continue;
1103         // apply correction for material of the current layer
1104         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1105         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1106         vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1107         vtrack->SetClIndex(ilayer,-1);
1108         modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1109         if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1110           vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1111         }
1112         if(constrain && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1113           fUseImproveKalman ?
1114             vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1115             vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1116         }
1117         ntracks[ilayer]++;
1118         continue;
1119       }
1120
1121       // track outside layer acceptance in z
1122       if (idet<0) continue;
1123
1124       //propagate to the intersection with the detector plane
1125       const AliITSdetector &det=layer.GetDetector(idet);
1126       new(&currenttrack2)  AliITStrackMI(currenttrack1);
1127       if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1128       if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1129       currenttrack1.SetDetectorIndex(idet);
1130       currenttrack2.SetDetectorIndex(idet);
1131       if(!LocalModuleCoord(ilayer,idet,&currenttrack1,xloc,zloc)) continue; // local module coords
1132
1133       //***************
1134       // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1135       //
1136       // road in global (rphi,z) [i.e. in tracking ref. system]
1137       Double_t zmin,zmax,ymin,ymax;
1138       if (!ComputeRoad(&currenttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1139
1140       // select clusters in road
1141       layer.SelectClusters(zmin,zmax,ymin,ymax);
1142       //********************
1143
1144       // Define criteria for track-cluster association
1145       Double_t msz = currenttrack1.GetSigmaZ2() +
1146         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1147         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1148         AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1149       Double_t msy = currenttrack1.GetSigmaY2() +
1150         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1151         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1152         AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1153       if (constrain) {
1154         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1155         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1156       }  else {
1157         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1158         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1159       }
1160       msz = 1./msz; // 1/RoadZ^2
1161       msy = 1./msy; // 1/RoadY^2
1162
1163       //
1164       //
1165       // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1166       //
1167       const AliITSRecPoint *cl=0;
1168       Int_t clidx=-1;
1169       Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1170       Bool_t deadzoneSPD=kFALSE;
1171       currenttrack = &currenttrack1;
1172
1173       // check if the road contains a dead zone
1174       Bool_t noClusters = kFALSE;
1175       if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1176       if (noClusters) AliDebug(2,"no clusters in road");
1177       Double_t dz=0.5*(zmax-zmin);
1178       Double_t dy=0.5*(ymax-ymin);
1179       Int_t dead = CheckDeadZone(&currenttrack1,ilayer,idet,dz,dy,noClusters);
1180       if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1181       // create a prolongation without clusters (check also if there are no clusters in the road)
1182       if (dead ||
1183           (noClusters &&
1184            AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1185         AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1186         updatetrack->SetClIndex(ilayer,-1);
1187         if (dead==0) {
1188           modstatus = 5; // no cls in road
1189         } else if (dead==1) {
1190           modstatus = 7; // holes in z in SPD
1191         } else if (dead==2 || dead==3 || dead==4) {
1192           modstatus = 2; // dead from OCDB
1193         }
1194         updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1195         // apply correction for material of the current layer
1196         CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1197         if (constrain) { // apply vertex constrain
1198           updatetrack->SetConstrain(constrain);
1199           Bool_t isPrim = kTRUE;
1200           if (ilayer<4) { // check that it's close to the vertex
1201             updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1202             if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1203                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1204                 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1205                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1206           }
1207           if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1208             fUseImproveKalman ?
1209               updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1210               updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1211           }
1212         }
1213         updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1214         if (dead) {
1215           if (dead==1) { // dead zone at z=0,+-7cm in SPD
1216             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1217             deadzoneSPD=kTRUE;
1218           } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1219             updatetrack->SetDeadZoneProbability(ilayer,1.);
1220           } else if (dead==4) { // at least a single dead channel from OCDB
1221             updatetrack->SetDeadZoneProbability(ilayer,0.);
1222           }
1223         }
1224         ntracks[ilayer]++;
1225       }
1226
1227       clidx=-1;
1228       // loop over clusters in the road
1229       while ((cl=layer.GetNextCluster(clidx))!=0) {
1230         if (ntracks[ilayer]>int(0.95*kMaxTr)) break; //space for skipped clusters
1231         Bool_t changedet =kFALSE;
1232         if (TMath::Abs(cl->GetQ())<1.e-13 && deadzoneSPD==kTRUE) continue;
1233         Int_t idetc=cl->GetDetectorIndex();
1234
1235         if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1236           // take into account misalignment (bring track to real detector plane)
1237           Double_t xTrOrig = currenttrack->GetX();
1238           if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1239           // a first cut on track-cluster distance
1240           if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1241                (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1242             {  // cluster not associated to track
1243               AliDebug(2,"not associated");
1244               // MvL: added here as well
1245               // bring track back to ideal detector plane
1246               currenttrack->Propagate(xTrOrig);
1247               continue;
1248             }
1249           // bring track back to ideal detector plane
1250           if (!currenttrack->Propagate(xTrOrig)) continue;
1251         } else {                                      // have to move track to cluster's detector
1252           const AliITSdetector &detc=layer.GetDetector(idetc);
1253           // a first cut on track-cluster distance
1254           Double_t y;
1255           if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1256           if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1257                (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1258             continue; // cluster not associated to track
1259           //
1260           new (&backuptrack) AliITStrackMI(currenttrack2);
1261           changedet = kTRUE;
1262           currenttrack =&currenttrack2;
1263           if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1264             new (currenttrack) AliITStrackMI(backuptrack);
1265             changedet = kFALSE;
1266             continue;
1267           }
1268           currenttrack->SetDetectorIndex(idetc);
1269           // Get again the budget to the primary vertex
1270           // for the current track being prolonged, if had to change detector
1271           //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1272         }
1273
1274         // calculate track-clusters chi2
1275         chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1276         // chi2 cut
1277         AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1278         if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1279           if (TMath::Abs(cl->GetQ())<1.e-13) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1280           if (ntracks[ilayer]>=kMaxTr) continue;
1281           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1282           updatetrack->SetClIndex(ilayer,-1);
1283           if (changedet) new (&currenttrack2) AliITStrackMI(backuptrack);
1284
1285           if (TMath::Abs(cl->GetQ())>1.e-13) { // real cluster
1286             if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1287               AliDebug(2,"update failed");
1288               continue;
1289             }
1290             updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1291             modstatus = 1; // found
1292           } else {             // virtual cluster in dead zone
1293             updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1294             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1295             modstatus = 7; // holes in z in SPD
1296           }
1297
1298           if (changedet) {
1299             Float_t xlocnewdet,zlocnewdet;
1300             if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1301               updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1302             }
1303           } else {
1304             updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1305           }
1306           if (cl->IsUsed()) updatetrack->IncrementNUsed();
1307
1308           // apply correction for material of the current layer
1309           CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1310
1311           if (constrain) { // apply vertex constrain
1312             updatetrack->SetConstrain(constrain);
1313             Bool_t isPrim = kTRUE;
1314             if (ilayer<4) { // check that it's close to the vertex
1315               updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1316               if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1317                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1318                   TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1319                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1320             }
1321             if (isPrim && AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1322               fUseImproveKalman ?
1323                 updatetrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1324                 updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1325             }
1326           } //apply vertex constrain
1327           ntracks[ilayer]++;
1328         }  // create new hypothesis
1329         else {
1330           AliDebug(2,"chi2 too large");
1331         }
1332
1333       } // loop over possible prolongations
1334
1335       // allow one prolongation without clusters
1336       if (constrain && itrack<=1 && TMath::Abs(currenttrack1.GetNSkipped())<1.e-13 && deadzoneSPD==kFALSE && ntracks[ilayer]<kMaxTr) {
1337         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1338         // apply correction for material of the current layer
1339         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1340         vtrack->SetClIndex(ilayer,-1);
1341         modstatus = 3; // skipped
1342         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1343         if(AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
1344           fUseImproveKalman ?
1345             vtrack->ImproveKalman(xyzVtx,ersVtx,xMSLrs,x2X0MSLrs,nMSLrs) :
1346             vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1347         }
1348         vtrack->IncrementNSkipped();
1349         ntracks[ilayer]++;
1350       }
1351
1352
1353     } // loop over tracks in layer ilayer+1
1354
1355     //loop over track candidates for the current layer
1356     //
1357     //
1358     Int_t accepted=0;
1359
1360     Int_t golden=0;
1361     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1362       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1363       if (normalizedchi2[itrack] <
1364           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1365       if (ilayer>4) {
1366         accepted++;
1367       } else {
1368         if (constrain) { // constrain
1369           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1370             accepted++;
1371         } else {        // !constrain
1372           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1373             accepted++;
1374         }
1375       }
1376     }
1377     // sort tracks by increasing normalized chi2
1378     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1379     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1380     if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1381     //    if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1382     if (ntracks[ilayer]>int(kMaxTr*0.9)) ntracks[ilayer]=int(kMaxTr*0.9);
1383   } // end loop over layers
1384
1385
1386   //
1387   // Now select tracks to be kept
1388   //
1389   Int_t max = constrain ? 20 : 5;
1390
1391   // tracks that reach layer 0 (SPD inner)
1392   for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1393     AliITStrackMI & track= tracks[0][nindexes[0][i]];
1394     if (track.GetNumberOfClusters()<2) continue;
1395     if (!constrain && track.GetNormChi2(0) >
1396         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1397       continue;
1398     }
1399     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1400   }
1401
1402   // tracks that reach layer 1 (SPD outer)
1403   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1404     AliITStrackMI & track= tracks[1][nindexes[1][i]];
1405     if (track.GetNumberOfClusters()<4) continue;
1406     if (!constrain && track.GetNormChi2(1) >
1407         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1408     if (constrain) track.IncrementNSkipped();
1409     if (!constrain) {
1410       track.SetD(0,track.GetD(GetX(),GetY()));
1411       track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1412       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1413         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1414       }
1415     }
1416     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1417   }
1418
1419   // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1420   if (!constrain){
1421     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1422       AliITStrackMI & track= tracks[2][nindexes[2][i]];
1423       if (track.GetNumberOfClusters()<3) continue;
1424       if (track.GetNormChi2(2) >
1425           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1426       track.SetD(0,track.GetD(GetX(),GetY()));
1427       track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1428       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1429         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1430       }
1431       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1432     }
1433   }
1434
1435   if (!constrain) {
1436     //
1437     // register best track of each layer - important for V0 finder
1438     //
1439     for (Int_t ilayer=0;ilayer<5;ilayer++){
1440       if (ntracks[ilayer]==0) continue;
1441       AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1442       if (track.GetNumberOfClusters()<1) continue;
1443       CookLabel(&track,0);
1444       bestarray->AddAt(new AliITStrackMI(track),ilayer);
1445     }
1446   }
1447   //
1448   // update TPC V0 information
1449   //
1450   if (otrack->GetESDtrack()->GetV0Index(0)>0){
1451     Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1452     for (Int_t i=0;i<3;i++){
1453       Int_t  index = otrack->GetESDtrack()->GetV0Index(i);
1454       if (index==0) break;
1455       AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1456       if (vertex->GetStatus()<0) continue;     // rejected V0
1457       //
1458       if (otrack->GetSign()>0) {
1459         vertex->SetIndex(0,esdindex);
1460       }
1461       else{
1462         vertex->SetIndex(1,esdindex);
1463       }
1464       //find nearest layer with track info
1465       Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
1466       Int_t nearestold  = GetNearestLayer(xrp);               //I.B.
1467       Int_t nearest     = nearestold;
1468       for (Int_t ilayer =nearest;ilayer<7;ilayer++){
1469         if (ntracks[nearest]==0){
1470           nearest = ilayer;
1471         }
1472       }
1473       //
1474       AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1475       if (nearestold<5&&nearest<5){
1476         Bool_t accept = track.GetNormChi2(nearest)<10;
1477         if (accept){
1478           if (track.GetSign()>0) {
1479             vertex->SetParamP(track);
1480             vertex->Update(fprimvertex);
1481             //vertex->SetIndex(0,track.fESDtrack->GetID());
1482             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1483           }else{
1484             vertex->SetParamN(track);
1485             vertex->Update(fprimvertex);
1486             //vertex->SetIndex(1,track.fESDtrack->GetID());
1487             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1488           }
1489           vertex->SetStatus(vertex->GetStatus()+1);
1490         }else{
1491           //vertex->SetStatus(-2);  // reject V0  - not enough clusters
1492         }
1493       }
1494     }
1495   }
1496
1497 }
1498 //------------------------------------------------------------------------
1499 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1500 {
1501   //--------------------------------------------------------------------
1502   //
1503   //
1504   return fgLayers[layer];
1505 }
1506 //------------------------------------------------------------------------
1507 AliITStrackerMI::AliITSlayer::AliITSlayer():
1508 fR(0),
1509 fPhiOffset(0),
1510 fNladders(0),
1511 fZOffset(0),
1512 fNdetectors(0),
1513 fDetectors(0),
1514 fN(0),
1515 fDy5(0),
1516 fDy10(0),
1517 fDy20(0),
1518 fClustersCs(0),
1519 fClusterIndexCs(0),
1520 fYcs(0),
1521 fZcs(0),
1522 fNcs(0),
1523 fCurrentSlice(-1),
1524 fZmin(0),
1525 fZmax(0),
1526 fYmin(0),
1527 fYmax(0),
1528 fI(0),
1529 fImax(0),
1530 fSkip(0),
1531 fAccepted(0),
1532 fRoad(0),
1533 fMaxSigmaClY(0),
1534 fMaxSigmaClZ(0),
1535 fNMaxSigmaCl(3)
1536 {
1537   //--------------------------------------------------------------------
1538   //default AliITSlayer constructor
1539   //--------------------------------------------------------------------
1540   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1541     fClusterWeight[i]=0;
1542     fClusterTracks[0][i]=-1;
1543     fClusterTracks[1][i]=-1;
1544     fClusterTracks[2][i]=-1;
1545     fClusterTracks[3][i]=-1;
1546     fY[i]=0;
1547     fZ[i]=0;
1548   }
1549   fYB[0]=0;
1550   fYB[1]=0;
1551
1552   for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1553     for (Int_t j1=0; j1<6; j1++) {
1554       fClusters5[j1][j]=0;
1555       fClusterIndex5[j1][j]=-1;
1556       fY5[j1][j]=0;
1557       fZ5[j1][j]=0;
1558       fN5[j1]=0;
1559       fBy5[j1][0]=0;
1560       fBy5[j1][1]=0;
1561     }
1562   }
1563
1564   for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1565     for (Int_t j1=0; j1<11; j1++) {
1566       fClusters10[j1][j]=0;
1567       fClusterIndex10[j1][j]=-1;
1568       fY10[j1][j]=0;
1569       fZ10[j1][j]=0;
1570       fN10[j1]=0;
1571       fBy10[j1][0]=0;
1572       fBy10[j1][1]=0;
1573     }
1574   }
1575
1576   for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1577     for (Int_t j1=0; j1<21; j1++) {
1578       fClusters20[j1][j]=0;
1579       fClusterIndex20[j1][j]=-1;
1580       fY20[j1][j]=0;
1581       fZ20[j1][j]=0;
1582       fN20[j1]=0;
1583       fBy20[j1][0]=0;
1584       fBy20[j1][1]=0;
1585     }
1586   }
1587   for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1588     fClusters[i]=NULL;
1589     fClusterIndex[i]=0;
1590   }
1591 }
1592 //------------------------------------------------------------------------
1593 AliITStrackerMI::AliITSlayer::
1594 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1595 fR(r),
1596 fPhiOffset(p),
1597 fNladders(nl),
1598 fZOffset(z),
1599 fNdetectors(nd),
1600 fDetectors(0),
1601 fN(0),
1602 fDy5(0),
1603 fDy10(0),
1604 fDy20(0),
1605 fClustersCs(0),
1606 fClusterIndexCs(0),
1607 fYcs(0),
1608 fZcs(0),
1609 fNcs(0),
1610 fCurrentSlice(-1),
1611 fZmin(0),
1612 fZmax(0),
1613 fYmin(0),
1614 fYmax(0),
1615 fI(0),
1616 fImax(0),
1617 fSkip(0),
1618 fAccepted(0),
1619 fRoad(0),
1620 fMaxSigmaClY(0),
1621 fMaxSigmaClZ(0),
1622 fNMaxSigmaCl(3) {
1623   //--------------------------------------------------------------------
1624   //main AliITSlayer constructor
1625   //--------------------------------------------------------------------
1626   fDetectors=new AliITSdetector[fNladders*fNdetectors];
1627   fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1628
1629   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1630     fClusterWeight[i]=0;
1631     fClusterTracks[0][i]=-1;
1632     fClusterTracks[1][i]=-1;
1633     fClusterTracks[2][i]=-1;
1634     fClusterTracks[3][i]=-1;
1635     fY[i]=0;
1636     fZ[i]=0;
1637   }
1638
1639   fYB[0]=0;
1640   fYB[1]=0;
1641
1642  for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer5; j++) {
1643     for (Int_t j1=0; j1<6; j1++) {
1644       fClusters5[j1][j]=0;
1645       fClusterIndex5[j1][j]=-1;
1646       fY5[j1][j]=0;
1647       fZ5[j1][j]=0;
1648       fN5[j1]=0;
1649       fBy5[j1][0]=0;
1650       fBy5[j1][1]=0;
1651     }
1652   }
1653
1654   for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer10; j++) {
1655     for (Int_t j1=0; j1<11; j1++) {
1656       fClusters10[j1][j]=0;
1657       fClusterIndex10[j1][j]=-1;
1658       fY10[j1][j]=0;
1659       fZ10[j1][j]=0;
1660       fN10[j1]=0;
1661       fBy10[j1][0]=0;
1662       fBy10[j1][1]=0;
1663     }
1664   }
1665
1666   for (Int_t j=0; j<AliITSRecoParam::kMaxClusterPerLayer20; j++) {
1667     for (Int_t j1=0; j1<21; j1++) {
1668       fClusters20[j1][j]=0;
1669       fClusterIndex20[j1][j]=-1;
1670       fY20[j1][j]=0;
1671       fZ20[j1][j]=0;
1672       fN20[j1]=0;
1673       fBy20[j1][0]=0;
1674       fBy20[j1][1]=0;
1675     }
1676   }
1677   for(Int_t i=0;i<AliITSRecoParam::kMaxClusterPerLayer;i++){
1678     fClusters[i]=NULL;
1679     fClusterIndex[i]=0;
1680   }
1681 }
1682 /*
1683 //------------------------------------------------------------------------
1684 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1685 fR(layer.fR),
1686 fPhiOffset(layer.fPhiOffset),
1687 fNladders(layer.fNladders),
1688 fZOffset(layer.fZOffset),
1689 fNdetectors(layer.fNdetectors),
1690 fDetectors(layer.fDetectors),
1691 fN(layer.fN),
1692 fDy5(layer.fDy5),
1693 fDy10(layer.fDy10),
1694 fDy20(layer.fDy20),
1695 fClustersCs(layer.fClustersCs),
1696 fClusterIndexCs(layer.fClusterIndexCs),
1697 fYcs(layer.fYcs),
1698 fZcs(layer.fZcs),
1699 fNcs(layer.fNcs),
1700 fCurrentSlice(layer.fCurrentSlice),
1701 fZmin(layer.fZmin),
1702 fZmax(layer.fZmax),
1703 fYmin(layer.fYmin),
1704 fYmax(layer.fYmax),
1705 fI(layer.fI),
1706 fImax(layer.fImax),
1707 fSkip(layer.fSkip),
1708 fAccepted(layer.fAccepted),
1709 fRoad(layer.fRoad),
1710 fMaxSigmaClY(layer.fMaxSigmaClY),
1711 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1712 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1713 {
1714   //Copy constructor
1715 }
1716 */
1717 //------------------------------------------------------------------------
1718 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1719   //--------------------------------------------------------------------
1720   // AliITSlayer destructor
1721   //--------------------------------------------------------------------
1722   delete [] fDetectors;
1723   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1724   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1725     fClusterWeight[i]=0;
1726     fClusterTracks[0][i]=-1;
1727     fClusterTracks[1][i]=-1;
1728     fClusterTracks[2][i]=-1;
1729     fClusterTracks[3][i]=-1;
1730   }
1731 }
1732 //------------------------------------------------------------------------
1733 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1734   //--------------------------------------------------------------------
1735   // This function removes loaded clusters
1736   //--------------------------------------------------------------------
1737   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1738   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1739     fClusterWeight[i]=0;
1740     fClusterTracks[0][i]=-1;
1741     fClusterTracks[1][i]=-1;
1742     fClusterTracks[2][i]=-1;
1743     fClusterTracks[3][i]=-1;
1744   }
1745
1746   fN=0;
1747   fI=0;
1748 }
1749 //------------------------------------------------------------------------
1750 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1751   //--------------------------------------------------------------------
1752   // This function reset weights of the clusters
1753   //--------------------------------------------------------------------
1754   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1755     fClusterWeight[i]=0;
1756     fClusterTracks[0][i]=-1;
1757     fClusterTracks[1][i]=-1;
1758     fClusterTracks[2][i]=-1;
1759     fClusterTracks[3][i]=-1;
1760   }
1761   for (Int_t i=0; i<fN;i++) {
1762     AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1763     if (cl&&cl->IsUsed()) cl->Use();
1764   }
1765
1766 }
1767 //------------------------------------------------------------------------
1768 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1769   //--------------------------------------------------------------------
1770   // This function calculates the road defined by the cluster density
1771   //--------------------------------------------------------------------
1772   Int_t n=0;
1773   for (Int_t i=0; i<fN; i++) {
1774      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1775   }
1776   if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1777 }
1778 //------------------------------------------------------------------------
1779 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1780   //--------------------------------------------------------------------
1781   //This function adds a cluster to this layer
1782   //--------------------------------------------------------------------
1783   if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1784     return 1;
1785   }
1786   fCurrentSlice=-1;
1787   fClusters[fN]=cl;
1788   fN++;
1789   AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1790   //AD
1791   Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1792   Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1793   if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1794   if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1795   if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1796   if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1797   //AD
1798   /*
1799   if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1800   if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1801   if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1802   if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1803   */
1804   return 0;
1805 }
1806 //------------------------------------------------------------------------
1807 void  AliITStrackerMI::AliITSlayer::SortClusters()
1808 {
1809   //
1810   //sort clusters
1811   //
1812   AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1813   Float_t *z                = new Float_t[fN];
1814   Int_t   * index           = new Int_t[fN];
1815   //
1816   fMaxSigmaClY=0.; //AD
1817   fMaxSigmaClZ=0.; //AD
1818
1819   for (Int_t i=0;i<fN;i++){
1820     z[i] = fClusters[i]->GetZ();
1821     // save largest errors in y and z for this layer
1822     fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1823     fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1824   }
1825   TMath::Sort(fN,z,index,kFALSE);
1826   for (Int_t i=0;i<fN;i++){
1827     clusters[i] = fClusters[index[i]];
1828   }
1829   //
1830   for (Int_t i=0;i<fN;i++){
1831     fClusters[i] = clusters[i];
1832     fZ[i]        = fClusters[i]->GetZ();
1833     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1834     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1835     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1836     fY[i] = y;
1837   }
1838   delete[] index;
1839   delete[] z;
1840   delete[] clusters;
1841   //
1842
1843   fYB[0]=10000000;
1844   fYB[1]=-10000000;
1845   for (Int_t i=0;i<fN;i++){
1846     if (fY[i]<fYB[0]) fYB[0]=fY[i];
1847     if (fY[i]>fYB[1]) fYB[1]=fY[i];
1848     fClusterIndex[i] = i;
1849   }
1850   //
1851   // fill slices
1852   fDy5 = (fYB[1]-fYB[0])/5.;
1853   fDy10 = (fYB[1]-fYB[0])/10.;
1854   fDy20 = (fYB[1]-fYB[0])/20.;
1855   for (Int_t i=0;i<6;i++)  fN5[i] =0;
1856   for (Int_t i=0;i<11;i++) fN10[i]=0;
1857   for (Int_t i=0;i<21;i++) fN20[i]=0;
1858   //
1859   for (Int_t i=0;i<6;i++) {fBy5[i][0] =  fYB[0]+(i-0.75)*fDy5; fBy5[i][1] =  fYB[0]+(i+0.75)*fDy5;}
1860   for (Int_t i=0;i<11;i++) {fBy10[i][0] =  fYB[0]+(i-0.75)*fDy10; fBy10[i][1] =  fYB[0]+(i+0.75)*fDy10;}
1861   for (Int_t i=0;i<21;i++) {fBy20[i][0] =  fYB[0]+(i-0.75)*fDy20; fBy20[i][1] =  fYB[0]+(i+0.75)*fDy20;}
1862   //
1863   //
1864   for (Int_t i=0;i<fN;i++)
1865     for (Int_t irot=-1;irot<=1;irot++){
1866       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1867       // slice 5
1868       for (Int_t slice=0; slice<6;slice++){
1869         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1870           fClusters5[slice][fN5[slice]] = fClusters[i];
1871           fY5[slice][fN5[slice]] = curY;
1872           fZ5[slice][fN5[slice]] = fZ[i];
1873           fClusterIndex5[slice][fN5[slice]]=i;
1874           fN5[slice]++;
1875         }
1876       }
1877       // slice 10
1878       for (Int_t slice=0; slice<11;slice++){
1879         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1880           fClusters10[slice][fN10[slice]] = fClusters[i];
1881           fY10[slice][fN10[slice]] = curY;
1882           fZ10[slice][fN10[slice]] = fZ[i];
1883           fClusterIndex10[slice][fN10[slice]]=i;
1884           fN10[slice]++;
1885         }
1886       }
1887       // slice 20
1888       for (Int_t slice=0; slice<21;slice++){
1889         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1890           fClusters20[slice][fN20[slice]] = fClusters[i];
1891           fY20[slice][fN20[slice]] = curY;
1892           fZ20[slice][fN20[slice]] = fZ[i];
1893           fClusterIndex20[slice][fN20[slice]]=i;
1894           fN20[slice]++;
1895         }
1896       }
1897     }
1898
1899   //
1900   // consistency check
1901   //
1902   for (Int_t i=0;i<fN-1;i++){
1903     if (fZ[i]>fZ[i+1]){
1904       printf("Bug\n");
1905     }
1906   }
1907   //
1908   for (Int_t slice=0;slice<21;slice++)
1909   for (Int_t i=0;i<fN20[slice]-1;i++){
1910     if (fZ20[slice][i]>fZ20[slice][i+1]){
1911       printf("Bug\n");
1912     }
1913   }
1914
1915
1916 }
1917 //------------------------------------------------------------------------
1918 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1919   //--------------------------------------------------------------------
1920   // This function returns the index of the nearest cluster
1921   //--------------------------------------------------------------------
1922   Int_t ncl=0;
1923   const Float_t *zcl;
1924   if (fCurrentSlice<0) {
1925     ncl = fN;
1926     zcl   = fZ;
1927   }
1928   else{
1929     ncl   = fNcs;
1930     zcl   = fZcs;;
1931   }
1932
1933   if (ncl==0) return 0;
1934   Int_t b=0, e=ncl-1, m=(b+e)/2;
1935   for (; b<e; m=(b+e)/2) {
1936     //    if (z > fClusters[m]->GetZ()) b=m+1;
1937     if (z > zcl[m]) b=m+1;
1938     else e=m;
1939   }
1940   return m;
1941 }
1942 //------------------------------------------------------------------------
1943 Bool_t AliITStrackerMI::ComputeRoad(AliITStrackMI* track,Int_t ilayer,Int_t idet,Double_t &zmin,Double_t &zmax,Double_t &ymin,Double_t &ymax) const {
1944   //--------------------------------------------------------------------
1945   // This function computes the rectangular road for this track
1946   //--------------------------------------------------------------------
1947
1948
1949   AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1950   // take into account the misalignment: propagate track to misaligned detector plane
1951   if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1952
1953   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1954                     TMath::Sqrt(track->GetSigmaZ2() +
1955                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1956                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1957                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1958   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1959                     TMath::Sqrt(track->GetSigmaY2() +
1960                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1961                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1962                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1963
1964   // track at boundary between detectors, enlarge road
1965   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1966   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1967        (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1968        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1969        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1970     Float_t tgl = TMath::Abs(track->GetTgl());
1971     if (tgl > 1.) tgl=1.;
1972     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1973     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1974     Float_t snp = TMath::Abs(track->GetSnp());
1975     if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1976     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1977   } // boundary
1978
1979   // add to the road a term (up to 2-3 mm) to deal with misalignments
1980   dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1981   dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1982
1983   Double_t r = fgLayers[ilayer].GetR();
1984   zmin = track->GetZ() - dz;
1985   zmax = track->GetZ() + dz;
1986   ymin = track->GetY() + r*det.GetPhi() - dy;
1987   ymax = track->GetY() + r*det.GetPhi() + dy;
1988
1989   // bring track back to idead detector plane
1990   if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1991
1992   return kTRUE;
1993 }
1994 //------------------------------------------------------------------------
1995 void AliITStrackerMI::AliITSlayer::
1996 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1997   //--------------------------------------------------------------------
1998   // This function sets the "window"
1999   //--------------------------------------------------------------------
2000
2001   Double_t circle=2*TMath::Pi()*fR;
2002   fYmin = ymin;
2003   fYmax = ymax;
2004   fZmin = zmin;
2005   fZmax = zmax;
2006   // AD
2007   // enlarge road in y by maximum cluster error on this layer (3 sigma)
2008   fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
2009   fYmax += fNMaxSigmaCl*fMaxSigmaClY;
2010   fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
2011   fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
2012
2013   Float_t ymiddle = (fYmin+fYmax)*0.5;
2014   if (ymiddle<fYB[0]) {
2015     fYmin+=circle; fYmax+=circle; ymiddle+=circle;
2016   } else if (ymiddle>fYB[1]) {
2017     fYmin-=circle; fYmax-=circle; ymiddle-=circle;
2018   }
2019
2020   //
2021   fCurrentSlice =-1;
2022   // defualt take all
2023   fClustersCs = fClusters;
2024   fClusterIndexCs = fClusterIndex;
2025   fYcs  = fY;
2026   fZcs  = fZ;
2027   fNcs  = fN;
2028   //
2029   //is in 20 slice?
2030   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
2031     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
2032     if (slice<0) slice=0;
2033     if (slice>20) slice=20;
2034     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
2035     if (isOK) {
2036       fCurrentSlice=slice;
2037       fClustersCs = fClusters20[fCurrentSlice];
2038       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
2039       fYcs  = fY20[fCurrentSlice];
2040       fZcs  = fZ20[fCurrentSlice];
2041       fNcs  = fN20[fCurrentSlice];
2042     }
2043   }
2044   //
2045   //is in 10 slice?
2046   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
2047     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
2048     if (slice<0) slice=0;
2049     if (slice>10) slice=10;
2050     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
2051     if (isOK) {
2052       fCurrentSlice=slice;
2053       fClustersCs = fClusters10[fCurrentSlice];
2054       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
2055       fYcs  = fY10[fCurrentSlice];
2056       fZcs  = fZ10[fCurrentSlice];
2057       fNcs  = fN10[fCurrentSlice];
2058     }
2059   }
2060   //
2061   //is in 5 slice?
2062   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
2063     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
2064     if (slice<0) slice=0;
2065     if (slice>5) slice=5;
2066     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
2067     if (isOK) {
2068       fCurrentSlice=slice;
2069       fClustersCs = fClusters5[fCurrentSlice];
2070       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
2071       fYcs  = fY5[fCurrentSlice];
2072       fZcs  = fZ5[fCurrentSlice];
2073       fNcs  = fN5[fCurrentSlice];
2074     }
2075   }
2076   //
2077   fI        = FindClusterIndex(fZmin);
2078   fImax     = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
2079   fSkip     = 0;
2080   fAccepted = 0;
2081
2082   return;
2083 }
2084 //------------------------------------------------------------------------
2085 Int_t AliITStrackerMI::AliITSlayer::
2086 FindDetectorIndex(Double_t phi, Double_t z) const {
2087   //--------------------------------------------------------------------
2088   //This function finds the detector crossed by the track
2089   //--------------------------------------------------------------------
2090   Double_t dphi;
2091   if (fZOffset<0)            // old geometry
2092     dphi = -(phi-fPhiOffset);
2093   else                       // new geometry
2094     dphi = phi-fPhiOffset;
2095
2096
2097   if      (dphi <  0) dphi += 2*TMath::Pi();
2098   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
2099   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
2100   if (np>=fNladders) np-=fNladders;
2101   if (np<0)          np+=fNladders;
2102
2103
2104   Double_t dz=fZOffset-z;
2105   Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
2106   Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
2107   if (nz>=fNdetectors || nz<0) {
2108     //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
2109     return -1;
2110   }
2111
2112   // ad hoc correction for 3rd ladder of SDD inner layer,
2113   // which is reversed (rotated by pi around local y)
2114   // this correction is OK only from AliITSv11Hybrid onwards
2115   if (GetR()>12. && GetR()<20.) { // SDD inner
2116     if(np==2) { // 3rd ladder
2117       Double_t posMod252[3];
2118       AliITSgeomTGeo::GetTranslation(252,posMod252);
2119       // check the Z coordinate of Mod 252: if negative
2120       // (old SDD geometry in AliITSv11Hybrid)
2121       // the swap of numeration whould be applied
2122       if(posMod252[2]<0.){
2123         nz = (fNdetectors-1) - nz;
2124       }
2125     }
2126   }
2127   //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
2128
2129
2130   return np*fNdetectors + nz;
2131 }
2132 //------------------------------------------------------------------------
2133 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
2134 {
2135   //--------------------------------------------------------------------
2136   // This function returns clusters within the "window"
2137   //--------------------------------------------------------------------
2138
2139   if (fCurrentSlice<0) {
2140     Double_t rpi2 = 2.*fR*TMath::Pi();
2141     for (Int_t i=fI; i<fImax; i++) {
2142       Double_t y = fY[i];
2143       Double_t z = fZ[i];
2144       if (fYmax<y) y -= rpi2;
2145       if (fYmin>y) y += rpi2;
2146       if (y<fYmin) continue;
2147       if (y>fYmax) continue;
2148       // AD
2149       // skip clusters that are in "extended" road but they
2150       // 3sigma error does not touch the original road
2151       if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
2152       if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
2153       //
2154       if (TMath::Abs(fClusters[i]->GetQ())<1.e-13 && fSkip==2) continue;
2155       ci=i;
2156       if (!test) fI=i+1;
2157       return fClusters[i];
2158     }
2159   } else {
2160     for (Int_t i=fI; i<fImax; i++) {
2161       if (fYcs[i]<fYmin) continue;
2162       if (fYcs[i]>fYmax) continue;
2163       if (TMath::Abs(fClustersCs[i]->GetQ())<1.e-13 && fSkip==2) continue;
2164       ci=fClusterIndexCs[i];
2165       if (!test) fI=i+1;
2166       return fClustersCs[i];
2167     }
2168   }
2169   return 0;
2170 }
2171 //------------------------------------------------------------------------
2172 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
2173 const {
2174   //--------------------------------------------------------------------
2175   // This function returns the layer thickness at this point (units X0)
2176   //--------------------------------------------------------------------
2177   Double_t d=0.0085;
2178   x0=AliITSRecoParam::GetX0Air();
2179   if (43<fR&&fR<45) { //SSD2
2180      Double_t dd=0.0034;
2181      d=dd;
2182      if (TMath::Abs(y-0.00)>3.40) d+=dd;
2183      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2184      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2185      for (Int_t i=0; i<12; i++) {
2186        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
2187           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2188           d+=0.0034;
2189           break;
2190        }
2191        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
2192           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2193           d+=0.0034;
2194           break;
2195        }
2196        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2197        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2198      }
2199   } else
2200   if (37<fR&&fR<41) { //SSD1
2201      Double_t dd=0.0034;
2202      d=dd;
2203      if (TMath::Abs(y-0.00)>3.40) d+=dd;
2204      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2205      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2206      for (Int_t i=0; i<11; i++) {
2207        if (TMath::Abs(z-3.9*i)<0.15) {
2208           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2209           d+=dd;
2210           break;
2211        }
2212        if (TMath::Abs(z+3.9*i)<0.15) {
2213           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2214           d+=dd;
2215           break;
2216        }
2217        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2218        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2219      }
2220   } else
2221   if (13<fR&&fR<26) { //SDD
2222      Double_t dd=0.0033;
2223      d=dd;
2224      if (TMath::Abs(y-0.00)>3.30) d+=dd;
2225
2226      if (TMath::Abs(y-1.80)<0.55) {
2227         d+=0.016;
2228         for (Int_t j=0; j<20; j++) {
2229           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2230           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2231         }
2232      }
2233      if (TMath::Abs(y+1.80)<0.55) {
2234         d+=0.016;
2235         for (Int_t j=0; j<20; j++) {
2236           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2237           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2238         }
2239      }
2240
2241      for (Int_t i=0; i<4; i++) {
2242        if (TMath::Abs(z-7.3*i)<0.60) {
2243           d+=dd;
2244           if (TMath::Abs(y-0.00)>3.30) d+=dd;
2245           break;
2246        }
2247        if (TMath::Abs(z+7.3*i)<0.60) {
2248           d+=dd;
2249           if (TMath::Abs(y-0.00)>3.30) d+=dd;
2250           break;
2251        }
2252      }
2253   } else
2254   if (6<fR&&fR<8) {   //SPD2
2255      Double_t dd=0.0063; x0=21.5;
2256      d=dd;
2257      if (TMath::Abs(y-3.08)>0.5) d+=dd;
2258      if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2259   } else
2260   if (3<fR&&fR<5) {   //SPD1
2261      Double_t dd=0.0063; x0=21.5;
2262      d=dd;
2263      if (TMath::Abs(y+0.21)>0.6) d+=dd;
2264      if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2265   }
2266
2267   return d;
2268 }
2269 //------------------------------------------------------------------------
2270 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2271 fR(det.fR),
2272 fRmisal(det.fRmisal),
2273 fPhi(det.fPhi),
2274 fSinPhi(det.fSinPhi),
2275 fCosPhi(det.fCosPhi),
2276 fYmin(det.fYmin),
2277 fYmax(det.fYmax),
2278 fZmin(det.fZmin),
2279 fZmax(det.fZmax),
2280 fIsBad(det.fIsBad),
2281 fNChips(det.fNChips),
2282 fChipIsBad(det.fChipIsBad)
2283 {
2284   //Copy constructor
2285 }
2286 //------------------------------------------------------------------------
2287 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2288                                                const AliITSDetTypeRec *detTypeRec)
2289 {
2290   //--------------------------------------------------------------------
2291   // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2292   //--------------------------------------------------------------------
2293
2294   // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2295   // while in the tracker they start from 0 for each layer
2296   for(Int_t il=0; il<ilayer; il++)
2297     idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2298
2299   Int_t detType;
2300   if (ilayer==0 || ilayer==1) {        // ----------  SPD
2301     detType = 0;
2302   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
2303     detType = 1;
2304   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
2305     detType = 2;
2306   } else {
2307     printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2308     return;
2309   }
2310
2311   // Get calibration from AliITSDetTypeRec
2312   AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2313   calib->SetModuleIndex(idet);
2314   AliITSCalibration *calibSPDdead = 0;
2315   if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2316   if (calib->IsBad() ||
2317       (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2318     {
2319       SetBad();
2320       //      printf("lay %d bad %d\n",ilayer,idet);
2321     }
2322
2323   // Get segmentation from AliITSDetTypeRec
2324   AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2325
2326   // Read info about bad chips
2327   fNChips = segm->GetMaximumChipIndex()+1;
2328   //printf("ilayer %d  detType %d idet %d fNChips %d %d  GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2329   if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2330   fChipIsBad = new Bool_t[fNChips];
2331   for (Int_t iCh=0;iCh<fNChips;iCh++) {
2332     fChipIsBad[iCh] = calib->IsChipBad(iCh);
2333     if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2334     //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2335   }
2336
2337   return;
2338 }
2339 //------------------------------------------------------------------------
2340 Double_t AliITStrackerMI::GetEffectiveThickness()
2341 {
2342   //--------------------------------------------------------------------
2343   // Returns the thickness between the current layer and the vertex (units X0)
2344   //--------------------------------------------------------------------
2345
2346   if(fUseTGeo!=0) {
2347     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2348     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2349     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2350   }
2351
2352   // beam pipe
2353   Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2354   Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2355
2356   // layers
2357   Double_t x0=0;
2358   Double_t xn=fgLayers[fI].GetR();
2359   for (Int_t i=0; i<fI; i++) {
2360     Double_t xi=fgLayers[i].GetR();
2361     Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2362     d+=dLayer*xi*xi;
2363   }
2364
2365   // shields
2366   if (fI>1) {
2367     Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2368     d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2369   }
2370   if (fI>3) {
2371     Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2372     d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2373   }
2374   return d/(xn*xn);
2375 }
2376
2377 //------------------------------------------------------------------------
2378 Int_t AliITStrackerMI::GetEffectiveThicknessLbyL(Double_t* xMS, Double_t* x2x0MS)
2379 {
2380   //--------------------------------------------------------------------
2381   // Returns the array of layers between the current layer and the vertex
2382   //--------------------------------------------------------------------
2383   //
2384   if(fUseTGeo!=0) {
2385     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2386     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2387     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2388   }
2389
2390   int nl = 0;
2391   double x0 = 0;
2392   for (int il=fI;il--;) {
2393     //
2394     if (il==3) {
2395       x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2396       xMS[nl++]  = AliITSRecoParam::GetrInsideShield(1);
2397     }
2398     else if (il==1) {
2399       x2x0MS[nl] = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2400       xMS[nl++]  = AliITSRecoParam::GetrInsideShield(0);
2401     }
2402     //
2403     x2x0MS[nl] = (fUseTGeo==0 ? fgLayers[il].GetThickness(0,0,x0) : fxOverX0Layer[il]);
2404     xMS[nl++]  = fgLayers[il].GetR();
2405     //
2406   }
2407   //
2408   // beam pipe
2409   x2x0MS[nl]  = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2410   xMS[nl++]  = AliITSRecoParam::GetrPipe();
2411   //
2412   return nl;
2413 }
2414
2415
2416 //------------------------------------------------------------------------
2417 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2418   //-------------------------------------------------------------------
2419   // This function returns number of clusters within the "window"
2420   //--------------------------------------------------------------------
2421   Int_t ncl=0;
2422   for (Int_t i=fI; i<fN; i++) {
2423     const AliITSRecPoint *c=fClusters[i];
2424     if (c->GetZ() > fZmax) break;
2425     if (c->IsUsed()) continue;
2426     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2427     Double_t y=fR*det.GetPhi() + c->GetY();
2428
2429     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2430     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2431
2432     if (y<fYmin) continue;
2433     if (y>fYmax) continue;
2434     ncl++;
2435   }
2436   return ncl;
2437 }
2438 //------------------------------------------------------------------------
2439 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2440                                 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2441 {
2442   //--------------------------------------------------------------------
2443   // This function refits the track "track" at the position "x" using
2444   // the clusters from "clusters"
2445   // If "extra"==kTRUE,
2446   //    the clusters from overlapped modules get attached to "track"
2447   // If "planeff"==kTRUE,
2448   //    special approach for plane efficiency evaluation is applyed
2449   //--------------------------------------------------------------------
2450
2451   Int_t index[AliITSgeomTGeo::kNLayers];
2452   Int_t k;
2453   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2454   Int_t nc=clusters->GetNumberOfClusters();
2455   for (k=0; k<nc; k++) {
2456     Int_t idx=clusters->GetClusterIndex(k);
2457     Int_t ilayer=(idx&0xf0000000)>>28;
2458     index[ilayer]=idx;
2459   }
2460
2461   return RefitAt(xx,track,index,extra,planeeff); // call the method below
2462 }
2463 //------------------------------------------------------------------------
2464 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2465                                 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2466 {
2467   //--------------------------------------------------------------------
2468   // This function refits the track "track" at the position "x" using
2469   // the clusters from array
2470   // If "extra"==kTRUE,
2471   //    the clusters from overlapped modules get attached to "track"
2472   // If "planeff"==kTRUE,
2473   //    special approach for plane efficiency evaluation is applyed
2474   //--------------------------------------------------------------------
2475   Int_t index[AliITSgeomTGeo::kNLayers];
2476   Int_t k;
2477   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2478   //
2479   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2480     index[k]=clusters[k];
2481   }
2482
2483   // special for cosmics and TPC prolonged tracks:
2484   // propagate to the innermost of:
2485   // - innermost layer crossed by the track
2486   // - innermost layer where a cluster was associated to the track
2487   static AliITSRecoParam *repa = NULL;
2488   if(!repa){
2489     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
2490     if(!repa){
2491       repa = AliITSRecoParam::GetHighFluxParam();
2492       AliWarning("Using default AliITSRecoParam class");
2493     }
2494   }
2495   Int_t evsp=repa->GetEventSpecie();
2496   ULong_t trStatus=0;
2497   if(track->GetESDtrack()) trStatus=track->GetStatus();
2498   Int_t innermostlayer=0;
2499   if((evsp&AliRecoParam::kCosmic) || (trStatus&AliESDtrack::kTPCin))  {
2500     innermostlayer=5;
2501     Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2502     for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2503       if( (drphi < (fgLayers[innermostlayer].GetR()+1.)) ||
2504           index[innermostlayer] >= 0 ) break;
2505     }
2506
2507     AliDebug(2,Form(" drphi  %f  innermost %d",drphi,innermostlayer));
2508   }
2509
2510   Int_t modstatus=1; // found
2511   Float_t xloc,zloc;
2512   Int_t from, to, step;
2513   if (xx > track->GetX()) {
2514       from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2515       step=+1;
2516   } else {
2517       from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2518       step=-1;
2519   }
2520   TString dir = (step>0 ? "outward" : "inward");
2521
2522   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2523      AliITSlayer &layer=fgLayers[ilayer];
2524      Double_t r=layer.GetR();
2525
2526      if (step<0 && xx>r) break;
2527
2528      // material between SSD and SDD, SDD and SPD
2529      Double_t hI=ilayer-0.5*step;
2530      if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2531        if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2532      if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2533        if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2534
2535
2536      Double_t oldGlobXYZ[3];
2537      if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2538
2539      // continue if we are already beyond this layer
2540      Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2541      if(step>0 && oldGlobR > r) continue; // going outward
2542      if(step<0 && oldGlobR < r) continue; // going inward
2543
2544      Double_t phi,z;
2545      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2546
2547      Int_t idet=layer.FindDetectorIndex(phi,z);
2548
2549      // check if we allow a prolongation without point for large-eta tracks
2550      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2551      if (skip==2) {
2552        modstatus = 4; // out in z
2553        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2554          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2555        }
2556        // cross layer
2557        // apply correction for material of the current layer
2558        // add time if going outward
2559        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2560        continue;
2561      }
2562
2563      if (idet<0) return kFALSE;
2564
2565
2566      const AliITSdetector &det=layer.GetDetector(idet);
2567      // only for ITS-SA tracks refit
2568      if (ilayer>1 && fTrackingPhase.Contains("RefitInward") && !(track->GetStatus()&AliESDtrack::kTPCin)) track->SetCheckInvariant(kFALSE);
2569      //
2570      if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
2571
2572      track->SetDetectorIndex(idet);
2573      if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2574
2575      Double_t dz,zmin,zmax,dy,ymin,ymax;
2576
2577      const AliITSRecPoint *clAcc=0;
2578      Double_t maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2579
2580      Int_t idx=index[ilayer];
2581      if (idx>=0) { // cluster in this layer
2582        modstatus = 6; // no refit
2583        const AliITSRecPoint *cl=(AliITSRecPoint *)GetCluster(idx);
2584        if (cl) {
2585          if (idet != cl->GetDetectorIndex()) {
2586            idet=cl->GetDetectorIndex();
2587            const AliITSdetector &detc=layer.GetDetector(idet);
2588            if (!track->Propagate(detc.GetPhi(),detc.GetR())) return kFALSE;
2589            track->SetDetectorIndex(idet);
2590            if(!LocalModuleCoord(ilayer,idet,track,xloc,zloc)) return kFALSE; // local module coords
2591          }
2592          Int_t cllayer = (idx & 0xf0000000) >> 28;;
2593          Double_t chi2=GetPredictedChi2MI(track,cl,cllayer);
2594          if (chi2<maxchi2) {
2595            clAcc=cl;
2596            maxchi2=chi2;
2597            modstatus = 1; // found
2598          } else {
2599             return kFALSE; //
2600          }
2601        }
2602      } else { // no cluster in this layer
2603        if (skip==1) {
2604          modstatus = 3; // skipped
2605          // Plane Eff determination:
2606          if (planeeff && ilayer==AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()) {
2607            if (IsOKForPlaneEff(track,clusters,ilayer))  // only adequate track for plane eff. evaluation
2608               UseTrackForPlaneEff(track,ilayer);
2609          }
2610        } else {
2611          modstatus = 5; // no cls in road
2612          // check dead
2613          if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2614          dz = 0.5*(zmax-zmin);
2615          dy = 0.5*(ymax-ymin);
2616          Int_t dead = CheckDeadZone(track,ilayer,idet,dz,dy,kTRUE);
2617          if (dead==1) modstatus = 7; // holes in z in SPD
2618          if (dead==2 || dead==3 || dead==4) modstatus = 2; // dead from OCDB
2619        }
2620      }
2621
2622      if (clAcc) {
2623        if (!UpdateMI(track,clAcc,maxchi2,idx)) return kFALSE;
2624        track->SetSampledEdx(clAcc->GetQ(),ilayer-2);
2625      }
2626      track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2627
2628
2629      if (extra && clAcc) { // search for extra clusters in overlapped modules
2630        AliITStrackV2 tmp(*track);
2631        if (!ComputeRoad(track,ilayer,idet,zmin,zmax,ymin,ymax)) return kFALSE;
2632        layer.SelectClusters(zmin,zmax,ymin,ymax);
2633
2634        const AliITSRecPoint *clExtra=0; Int_t ci=-1,cci=-1;
2635        Int_t idetExtra=-1;
2636        maxchi2=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
2637        Double_t tolerance=0.1;
2638        while ((clExtra=layer.GetNextCluster(ci))!=0) {
2639          // only clusters in another module! (overlaps)
2640          idetExtra = clExtra->GetDetectorIndex();
2641          if (idet == idetExtra) continue;
2642
2643          const AliITSdetector &detx=layer.GetDetector(idetExtra);
2644
2645          if (!tmp.Propagate(detx.GetPhi(),detx.GetR()+clExtra->GetX())) continue;
2646          if (TMath::Abs(tmp.GetZ() - clExtra->GetZ()) > tolerance) continue;
2647          if (TMath::Abs(tmp.GetY() - clExtra->GetY()) > tolerance) continue;
2648          if (!tmp.Propagate(detx.GetPhi(),detx.GetR())) continue;
2649
2650          Double_t chi2=tmp.GetPredictedChi2(clExtra);
2651          if (chi2<maxchi2) { maxchi2=chi2; cci=ci; }
2652        }
2653        if (cci>=0) {
2654          track->SetExtraCluster(ilayer,(ilayer<<28)+cci);
2655          track->SetExtraModule(ilayer,idetExtra);
2656        }
2657      } // end search for extra clusters in overlapped modules
2658
2659      // Correct for material of the current layer
2660      // cross material
2661      // add time if going outward
2662      if(!CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir)) return kFALSE;
2663      track->SetCheckInvariant(kTRUE);
2664   } // end loop on layers
2665
2666   if (!track->PropagateTo(xx,0.,0.)) return kFALSE;
2667
2668   return kTRUE;
2669 }
2670 //------------------------------------------------------------------------
2671 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
2672 {
2673   //
2674   // calculate normalized chi2
2675   //  return NormalizedChi2(track,0);
2676   Float_t chi2 = 0;
2677   Float_t sum=0;
2678   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2679   //  track->fdEdxMismatch=0;
2680   Float_t dedxmismatch =0;
2681   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
2682   if (mode<100){
2683     for (Int_t i = 0;i<6;i++){
2684       if (track->GetClIndex(i)>=0){
2685         Float_t cerry, cerrz;
2686         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2687         else
2688           { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2689         cerry*=cerry;
2690         cerrz*=cerrz;
2691         Float_t cchi2 = (track->GetDy(i)*track->GetDy(i)/cerry)+(track->GetDz(i)*track->GetDz(i)/cerrz);
2692         if (i>1 && AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(i)) {
2693           Float_t ratio = track->GetNormQ(i)/track->GetExpQ();
2694           if (ratio<0.5) {
2695             cchi2+=(0.5-ratio)*10.;
2696             //track->fdEdxMismatch+=(0.5-ratio)*10.;
2697             dedxmismatch+=(0.5-ratio)*10.;
2698           }
2699         }
2700         if (i<2 ||i>3){
2701           AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster( track->GetClIndex(i));
2702           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
2703           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
2704           if (i<2) chi2+=2*cl->GetDeltaProbability();
2705         }
2706         chi2+=cchi2;
2707         sum++;
2708       }
2709     }
2710     if (TMath::Abs(track->GetdEdxMismatch()-dedxmismatch)>0.0001){
2711       track->SetdEdxMismatch(dedxmismatch);
2712     }
2713   }
2714   else{
2715     for (Int_t i = 0;i<4;i++){
2716       if (track->GetClIndex(i)>=0){
2717         Float_t cerry, cerrz;
2718         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2719         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2720         cerry*=cerry;
2721         cerrz*=cerrz;
2722         chi2+= (track->GetDy(i)*track->GetDy(i)/cerry);
2723         chi2+= (track->GetDz(i)*track->GetDz(i)/cerrz);
2724         sum++;
2725       }
2726     }
2727     for (Int_t i = 4;i<6;i++){
2728       if (track->GetClIndex(i)>=0){
2729         Float_t cerry, cerrz;
2730         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
2731         else { cerry= track->GetSigmaY(i); cerrz = track->GetSigmaZ(i);}
2732         cerry*=cerry;
2733         cerrz*=cerrz;
2734         Float_t cerryb, cerrzb;
2735         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
2736         else { cerryb= track->GetSigmaY(i+6); cerrzb = track->GetSigmaZ(i+6);}
2737         cerryb*=cerryb;
2738         cerrzb*=cerrzb;
2739         chi2+= TMath::Min((track->GetDy(i+6)*track->GetDy(i+6)/cerryb),track->GetDy(i)*track->GetDy(i)/cerry);
2740         chi2+= TMath::Min((track->GetDz(i+6)*track->GetDz(i+6)/cerrzb),track->GetDz(i)*track->GetDz(i)/cerrz);
2741         sum++;
2742       }
2743     }
2744   }
2745   if (track->GetESDtrack()->GetTPCsignal()>85){
2746     Float_t ratio = track->GetdEdx()/track->GetESDtrack()->GetTPCsignal();
2747     if (ratio<0.5) {
2748       chi2+=(0.5-ratio)*5.;
2749     }
2750     if (ratio>2){
2751       chi2+=(ratio-2.0)*3;
2752     }
2753   }
2754   //
2755   Double_t match = TMath::Sqrt(track->GetChi22());
2756   if (track->GetConstrain())  match/=track->GetNumberOfClusters();
2757   if (!track->GetConstrain()) {
2758     if (track->GetNumberOfClusters()>2) {
2759       match/=track->GetNumberOfClusters()-2.;
2760     } else {
2761       match=0;
2762     }
2763   }
2764   if (match<0) match=0;
2765
2766   // penalty factor for missing points (NDeadZone>0), but no penalty
2767   // for layer with deadZoneProb close to 1 (either we wanted to skip layer
2768   // or there is a dead from OCDB)
2769   Float_t deadzonefactor = 0.;
2770   if (track->GetNDeadZone()>0.) {
2771     Int_t sumDeadZoneProbability=0;
2772     for(Int_t ilay=0;ilay<6;ilay++) {
2773       if(track->GetDeadZoneProbability(ilay)>0.) sumDeadZoneProbability++;
2774     }
2775     Int_t nDeadZoneWithProbNot1=(Int_t)(track->GetNDeadZone())-sumDeadZoneProbability;
2776     if(nDeadZoneWithProbNot1>0) {
2777       Float_t deadZoneProbability = track->GetNDeadZone()-(Float_t)sumDeadZoneProbability;
2778       AliDebug(2,Form("nDeadZone %f sumDZProbability %d nDZWithProbNot1 %d deadZoneProb %f\n",track->GetNDeadZone(),sumDeadZoneProbability,nDeadZoneWithProbNot1,deadZoneProbability));
2779       deadZoneProbability /= (Float_t)nDeadZoneWithProbNot1;
2780       Float_t one = 1.;
2781       deadZoneProbability = TMath::Min(deadZoneProbability,one);
2782       deadzonefactor = 3.*(1.1-deadZoneProbability);
2783     }
2784   }
2785
2786   Double_t normchi2 = 2*track->GetNSkipped()+match+deadzonefactor+(1+(2*track->GetNSkipped()+deadzonefactor)/track->GetNumberOfClusters())*
2787     (chi2)/TMath::Max(double(sum-track->GetNSkipped()),
2788                                 1./(1.+track->GetNSkipped()));
2789   AliDebug(2,Form("match %f deadzonefactor %f chi2 %f sum %f skipped %f\n",match,deadzonefactor,chi2,sum,track->GetNSkipped()));
2790   AliDebug(2,Form("NormChi2 %f  cls %d\n",normchi2,track->GetNumberOfClusters()));
2791   return normchi2;
2792 }
2793 //------------------------------------------------------------------------
2794 Double_t AliITStrackerMI::GetMatchingChi2(const AliITStrackMI * track1,const AliITStrackMI * track2)
2795 {
2796   //
2797   // return matching chi2 between two tracks
2798   Double_t largeChi2=1000.;
2799
2800   AliITStrackMI track3(*track2);
2801   if (!track3.Propagate(track1->GetAlpha(),track1->GetX())) return largeChi2;
2802   TMatrixD vec(5,1);
2803   vec(0,0)=track1->GetY()   - track3.GetY();
2804   vec(1,0)=track1->GetZ()   - track3.GetZ();
2805   vec(2,0)=track1->GetSnp() - track3.GetSnp();
2806   vec(3,0)=track1->GetTgl() - track3.GetTgl();
2807   vec(4,0)=track1->GetSigned1Pt() - track3.GetSigned1Pt();
2808   //
2809   TMatrixD cov(5,5);
2810   cov(0,0) = track1->GetSigmaY2()+track3.GetSigmaY2();
2811   cov(1,1) = track1->GetSigmaZ2()+track3.GetSigmaZ2();
2812   cov(2,2) = track1->GetSigmaSnp2()+track3.GetSigmaSnp2();
2813   cov(3,3) = track1->GetSigmaTgl2()+track3.GetSigmaTgl2();
2814   cov(4,4) = track1->GetSigma1Pt2()+track3.GetSigma1Pt2();
2815
2816   cov(0,1)=cov(1,0) = track1->GetSigmaZY()+track3.GetSigmaZY();
2817   cov(0,2)=cov(2,0) = track1->GetSigmaSnpY()+track3.GetSigmaSnpY();
2818   cov(0,3)=cov(3,0) = track1->GetSigmaTglY()+track3.GetSigmaTglY();
2819   cov(0,4)=cov(4,0) = track1->GetSigma1PtY()+track3.GetSigma1PtY();
2820   //
2821   cov(1,2)=cov(2,1) = track1->GetSigmaSnpZ()+track3.GetSigmaSnpZ();
2822   cov(1,3)=cov(3,1) = track1->GetSigmaTglZ()+track3.GetSigmaTglZ();
2823   cov(1,4)=cov(4,1) = track1->GetSigma1PtZ()+track3.GetSigma1PtZ();
2824   //
2825   cov(2,3)=cov(3,2) = track1->GetSigmaTglSnp()+track3.GetSigmaTglSnp();
2826   cov(2,4)=cov(4,2) = track1->GetSigma1PtSnp()+track3.GetSigma1PtSnp();
2827   //
2828   cov(3,4)=cov(4,3) = track1->GetSigma1PtTgl()+track3.GetSigma1PtTgl();
2829
2830   cov.Invert();
2831   TMatrixD vec2(cov,TMatrixD::kMult,vec);
2832   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
2833   return chi2(0,0);
2834 }
2835 //------------------------------------------------------------------------
2836 Double_t  AliITStrackerMI::GetSPDDeadZoneProbability(Double_t zpos, Double_t zerr) const
2837 {
2838   //
2839   //  return probability that given point (characterized by z position and error)
2840   //  is in SPD dead zone
2841   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
2842   //
2843   Double_t probability = 0.;
2844   Double_t nearestz = 0.,distz=0.;
2845   Int_t    nearestzone = -1;
2846   Double_t mindistz = 1000.;
2847
2848   // find closest dead zone
2849   for (Int_t i=0; i<3; i++) {
2850     distz=TMath::Abs(zpos-0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]));
2851     if (distz<mindistz) {
2852       nearestzone=i;
2853       nearestz=0.5*(fSPDdetzcentre[i]+fSPDdetzcentre[i+1]);
2854       mindistz=distz;
2855     }
2856   }
2857
2858   // too far from dead zone
2859   if (TMath::Abs(zpos-nearestz)>0.25+3.*zerr) return probability;
2860
2861
2862   Double_t zmin, zmax;
2863   if (nearestzone==0) { // dead zone at z = -7
2864     zmin = fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2865     zmax = fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2866   } else if (nearestzone==1) { // dead zone at z = 0
2867     zmin = fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2868     zmax = fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2869   } else if (nearestzone==2) { // dead zone at z = +7
2870     zmin = fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength();
2871     zmax = fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength();
2872   } else {
2873     zmin = 0.;
2874     zmax = 0.;
2875   }
2876   // probability that the true z is in the range [zmin,zmax] (i.e. inside
2877   // dead zone)
2878   probability = 0.5*( AliMathBase::ErfFast((zpos-zmin)/zerr/TMath::Sqrt(2.)) -
2879                       AliMathBase::ErfFast((zpos-zmax)/zerr/TMath::Sqrt(2.)) );
2880   AliDebug(2,Form("zpos %f +- %f nearestzone %d  zmin zmax %f %f prob %f\n",zpos,zerr,nearestzone,zmin,zmax,probability));
2881   return probability;
2882 }
2883 //------------------------------------------------------------------------
2884 Double_t AliITStrackerMI::GetTruncatedChi2(const AliITStrackMI * track, Float_t fac)
2885 {
2886   //
2887   // calculate normalized chi2
2888   Float_t chi2[6];
2889   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
2890   Float_t ncl = 0;
2891   for (Int_t i = 0;i<6;i++){
2892     if (TMath::Abs(track->GetDy(i))>0){
2893       chi2[i]= (track->GetDy(i)/erry[i])*(track->GetDy(i)/erry[i]);
2894       chi2[i]+= (track->GetDz(i)/errz[i])*(track->GetDz(i)/errz[i]);
2895       ncl++;
2896     }
2897     else{chi2[i]=10000;}
2898   }
2899   Int_t index[6];
2900   TMath::Sort(6,chi2,index,kFALSE);
2901   Float_t max = float(ncl)*fac-1.;
2902   Float_t sumchi=0, sumweight=0;
2903   for (Int_t i=0;i<max+1;i++){
2904     Float_t weight = (i<max)?1.:(max+1.-i);
2905     sumchi+=weight*chi2[index[i]];
2906     sumweight+=weight;
2907   }
2908   Double_t normchi2 = sumchi/sumweight;
2909   return normchi2;
2910 }
2911 //------------------------------------------------------------------------
2912 Double_t AliITStrackerMI::GetInterpolatedChi2(const AliITStrackMI * forwardtrack,const AliITStrackMI * backtrack)
2913 {
2914   //
2915   // calculate normalized chi2
2916   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
2917   Int_t npoints = 0;
2918   Double_t res =0;
2919   for (Int_t i=0;i<6;i++){
2920     if ( (backtrack->GetSigmaY(i)<0.000000001) || (forwardtrack->GetSigmaY(i)<0.000000001)) continue;
2921     Double_t sy1 = forwardtrack->GetSigmaY(i);
2922     Double_t sz1 = forwardtrack->GetSigmaZ(i);
2923     Double_t sy2 = backtrack->GetSigmaY(i);
2924     Double_t sz2 = backtrack->GetSigmaZ(i);
2925     if (i<2){ sy2=1000.;sz2=1000;}
2926     //
2927     Double_t dy0 = (forwardtrack->GetDy(i)/(sy1*sy1) +backtrack->GetDy(i)/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
2928     Double_t dz0 = (forwardtrack->GetDz(i)/(sz1*sz1) +backtrack->GetDz(i)/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
2929     //
2930     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
2931     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
2932     //
2933     res+= nz0*nz0+ny0*ny0;
2934     npoints++;
2935   }
2936   if (npoints>1) return
2937                    TMath::Max(0.3*forwardtrack->OneOverPt()-0.5,0.)+
2938                    //2*forwardtrack->fNUsed+
2939                    res/TMath::Max(double(npoints-forwardtrack->GetNSkipped()),
2940                                   1./(1.+forwardtrack->GetNSkipped()));
2941   return 1000;
2942 }
2943 //------------------------------------------------------------------------
2944 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
2945   //--------------------------------------------------------------------
2946   //       Return pointer to a given cluster
2947   //--------------------------------------------------------------------
2948   Int_t l=(index & 0xf0000000) >> 28;
2949   Int_t c=(index & 0x0fffffff) >> 00;
2950   return fgLayers[l].GetWeight(c);
2951 }
2952 //------------------------------------------------------------------------
2953 void AliITStrackerMI::RegisterClusterTracks(const AliITStrackMI* track,Int_t id)
2954 {
2955   //---------------------------------------------
2956   // register track to the list
2957   //
2958   if (track->GetESDtrack()->GetKinkIndex(0)!=0) return;  //don't register kink tracks
2959   //
2960   //
2961   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2962     Int_t index = track->GetClusterIndex(icluster);
2963     Int_t l=(index & 0xf0000000) >> 28;
2964     Int_t c=(index & 0x0fffffff) >> 00;
2965     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2966     for (Int_t itrack=0;itrack<4;itrack++){
2967       if (fgLayers[l].GetClusterTracks(itrack,c)<0){
2968         fgLayers[l].SetClusterTracks(itrack,c,id);
2969         break;
2970       }
2971     }
2972   }
2973 }
2974 //------------------------------------------------------------------------
2975 void AliITStrackerMI::UnRegisterClusterTracks(const AliITStrackMI* track, Int_t id)
2976 {
2977   //---------------------------------------------
2978   // unregister track from the list
2979   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
2980     Int_t index = track->GetClusterIndex(icluster);
2981     Int_t l=(index & 0xf0000000) >> 28;
2982     Int_t c=(index & 0x0fffffff) >> 00;
2983     if (c>fgLayers[l].GetNumberOfClusters()) continue;
2984     for (Int_t itrack=0;itrack<4;itrack++){
2985       if (fgLayers[l].GetClusterTracks(itrack,c)==id){
2986         fgLayers[l].SetClusterTracks(itrack,c,-1);
2987       }
2988     }
2989   }
2990 }
2991 //------------------------------------------------------------------------
2992 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSRecPoint *clist[6])
2993 {
2994   //-------------------------------------------------------------
2995   //get number of shared clusters
2996   //-------------------------------------------------------------
2997   Float_t shared=0;
2998   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
2999   // mean  number of clusters
3000   Float_t *ny = GetNy(id), *nz = GetNz(id);
3001
3002
3003   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
3004     Int_t index = track->GetClusterIndex(icluster);
3005     Int_t l=(index & 0xf0000000) >> 28;
3006     Int_t c=(index & 0x0fffffff) >> 00;
3007     if (c>fgLayers[l].GetNumberOfClusters()) continue;
3008     // if (ny[l]<1.e-13){
3009     //   printf("problem\n");
3010     // }
3011     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3012     Float_t weight=1;
3013     //
3014     Float_t deltan = 0;
3015     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3016     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3017       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3018     if (l<2 || l>3){
3019       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3020     }
3021     else{
3022       deltan = (cl->GetNz()-nz[l]);
3023     }
3024     if (deltan>2.0) continue;  // extended - highly probable shared cluster
3025     weight = 2./TMath::Max(3.+deltan,2.);
3026     //
3027     for (Int_t itrack=0;itrack<4;itrack++){
3028       if (fgLayers[l].GetClusterTracks(itrack,c)>=0 && fgLayers[l].GetClusterTracks(itrack,c)!=id){
3029         list[l]=index;
3030         clist[l] = (AliITSRecPoint*)GetCluster(index);
3031         track->SetSharedWeight(l,weight);
3032         shared+=weight;
3033         break;
3034       }
3035     }
3036   }
3037   track->SetNUsed(shared);
3038   return shared;
3039 }
3040 //------------------------------------------------------------------------
3041 Int_t AliITStrackerMI::GetOverlapTrack(const AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
3042 {
3043   //
3044   // find first shared track
3045   //
3046   // mean  number of clusters
3047   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
3048   //
3049   for (Int_t i=0;i<6;i++) overlist[i]=-1;
3050   Int_t sharedtrack=100000;
3051   Int_t tracks[24],trackindex=0;
3052   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
3053   //
3054   for (Int_t icluster=0;icluster<6;icluster++){
3055     if (clusterlist[icluster]<0) continue;
3056     Int_t index = clusterlist[icluster];
3057     Int_t l=(index & 0xf0000000) >> 28;
3058     Int_t c=(index & 0x0fffffff) >> 00;
3059     // if (ny[l]<1.e-13){
3060     //   printf("problem\n");
3061     // }
3062     if (c>fgLayers[l].GetNumberOfClusters()) continue;
3063     //if (l>3) continue;
3064     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3065     //
3066     Float_t deltan = 0;
3067     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
3068     if (l>2&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(l))
3069       if (track->GetNormQ(l)/track->GetExpQ()>3.5) continue;
3070     if (l<2 || l>3){
3071       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
3072     }
3073     else{
3074       deltan = (cl->GetNz()-nz[l]);
3075     }
3076     if (deltan>2.0) continue;  // extended - highly probable shared cluster
3077     //
3078     for (Int_t itrack=3;itrack>=0;itrack--){
3079       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3080       if (fgLayers[l].GetClusterTracks(itrack,c)!=trackID){
3081        tracks[trackindex]  = fgLayers[l].GetClusterTracks(itrack,c);
3082        trackindex++;
3083       }
3084     }
3085   }
3086   if (trackindex==0) return -1;
3087   if (trackindex==1){
3088     sharedtrack = tracks[0];
3089   }else{
3090     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
3091     else{
3092       //
3093       Int_t tracks2[24], cluster[24];
3094       for (Int_t i=0;i<24;i++){ tracks2[i]=-1; cluster[i]=0;}
3095       Int_t index =0;
3096       //
3097       for (Int_t i=0;i<trackindex;i++){
3098         if (tracks[i]<0) continue;
3099         tracks2[index] = tracks[i];
3100         cluster[index]++;
3101         for (Int_t j=i+1;j<trackindex;j++){
3102           if (tracks[j]<0) continue;
3103           if (tracks[j]==tracks[i]){
3104             cluster[index]++;
3105             tracks[j]=-1;
3106           }
3107         }
3108         index++;
3109       }
3110       Int_t max=0;
3111       for (Int_t i=0;i<index;i++){
3112         if (cluster[index]>max) {
3113           sharedtrack=tracks2[index];
3114           max=cluster[index];
3115         }
3116       }
3117     }
3118   }
3119
3120   if (sharedtrack>=100000) return -1;
3121   //
3122   // make list of overlaps
3123   shared =0;
3124   for (Int_t icluster=0;icluster<6;icluster++){
3125     if (clusterlist[icluster]<0) continue;
3126     Int_t index = clusterlist[icluster];
3127     Int_t l=(index & 0xf0000000) >> 28;
3128     Int_t c=(index & 0x0fffffff) >> 00;
3129     if (c>fgLayers[l].GetNumberOfClusters()) continue;
3130     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(index);
3131     if (l==0 || l==1){
3132       if (cl->GetNy()>2) continue;
3133       if (cl->GetNz()>2) continue;
3134     }
3135     if (l==4 || l==5){
3136       if (cl->GetNy()>3) continue;
3137       if (cl->GetNz()>3) continue;
3138     }
3139     //
3140     for (Int_t itrack=3;itrack>=0;itrack--){
3141       if (fgLayers[l].GetClusterTracks(itrack,c)<0) continue;
3142       if (fgLayers[l].GetClusterTracks(itrack,c)==sharedtrack){
3143         overlist[l]=index;
3144         shared++;
3145       }
3146     }
3147   }
3148   return sharedtrack;
3149 }
3150 //------------------------------------------------------------------------
3151 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1,AliITStrackMI* original){
3152   //
3153   // try to find track hypothesys without conflicts
3154   // with minimal chi2;
3155   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
3156   Int_t entries1 = arr1->GetEntriesFast();
3157   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
3158   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
3159   Int_t entries2 = arr2->GetEntriesFast();
3160   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
3161   //
3162   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
3163   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
3164   if (track10->Pt()>0.5+track20->Pt()) return track10;
3165   //
3166   for (Int_t itrack=0;itrack<entries1;itrack++){
3167     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3168     UnRegisterClusterTracks(track,trackID1);
3169   }
3170   //
3171   for (Int_t itrack=0;itrack<entries2;itrack++){
3172     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3173     UnRegisterClusterTracks(track,trackID2);
3174   }
3175   Int_t index1=0;
3176   Int_t index2=0;
3177   Float_t maxconflicts=6;
3178   Double_t maxchi2 =1000.;
3179   //
3180   // get weight of hypothesys - using best hypothesy
3181   Double_t w1,w2;
3182
3183   Int_t list1[6],list2[6];
3184   AliITSRecPoint *clist1[6], *clist2[6] ;
3185   RegisterClusterTracks(track10,trackID1);
3186   RegisterClusterTracks(track20,trackID2);
3187   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
3188   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
3189   UnRegisterClusterTracks(track10,trackID1);
3190   UnRegisterClusterTracks(track20,trackID2);
3191   //
3192   // normalized chi2
3193   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
3194   Float_t nerry[6],nerrz[6];
3195   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
3196   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
3197   for (Int_t i=0;i<6;i++){
3198      if ( (erry1[i]>0) && (erry2[i]>0)) {
3199        nerry[i] = TMath::Min(erry1[i],erry2[i]);
3200        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
3201      }else{
3202        nerry[i] = TMath::Max(erry1[i],erry2[i]);
3203        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
3204      }
3205      if (TMath::Abs(track10->GetDy(i))>0.000000000000001){
3206        chi21 += track10->GetDy(i)*track10->GetDy(i)/(nerry[i]*nerry[i]);
3207        chi21 += track10->GetDz(i)*track10->GetDz(i)/(nerrz[i]*nerrz[i]);
3208        ncl1++;
3209      }
3210      if (TMath::Abs(track20->GetDy(i))>0.000000000000001){
3211        chi22 += track20->GetDy(i)*track20->GetDy(i)/(nerry[i]*nerry[i]);
3212        chi22 += track20->GetDz(i)*track20->GetDz(i)/(nerrz[i]*nerrz[i]);
3213        ncl2++;
3214      }
3215   }
3216   chi21/=ncl1;
3217   chi22/=ncl2;
3218   //
3219   //
3220   Float_t d1 = TMath::Sqrt(track10->GetD(0)*track10->GetD(0)+track10->GetD(1)*track10->GetD(1))+0.1;
3221   Float_t d2 = TMath::Sqrt(track20->GetD(0)*track20->GetD(0)+track20->GetD(1)*track20->GetD(1))+0.1;
3222   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
3223   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
3224   //
3225   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
3226         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
3227         +1.*track10->Pt()/(track10->Pt()+track20->Pt())
3228         );
3229   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
3230         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
3231         +1.*track20->Pt()/(track10->Pt()+track20->Pt())
3232         );
3233
3234   Double_t sumw = w1+w2;
3235   w1/=sumw;
3236   w2/=sumw;
3237   if (w1<w2*0.5) {
3238     w1 = (d2+0.5)/(d1+d2+1.);
3239     w2 = (d1+0.5)/(d1+d2+1.);
3240   }
3241   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
3242   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
3243   //
3244   // get pair of "best" hypothesys
3245   //
3246   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
3247   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
3248
3249   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3250     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
3251     //if (track1->fFakeRatio>0) continue;
3252     RegisterClusterTracks(track1,trackID1);
3253     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
3254       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
3255
3256       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
3257       //if (track2->fFakeRatio>0) continue;
3258       Float_t nskipped=0;
3259       RegisterClusterTracks(track2,trackID2);
3260       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
3261       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
3262       UnRegisterClusterTracks(track2,trackID2);
3263       //
3264       if (track1->GetConstrain()) nskipped+=w1*track1->GetNSkipped();
3265       if (track2->GetConstrain()) nskipped+=w2*track2->GetNSkipped();
3266       if (nskipped>0.5) continue;
3267       //
3268       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
3269       if (conflict1+1<cconflict1) continue;
3270       if (conflict2+1<cconflict2) continue;
3271       Float_t conflict=0;
3272       Float_t sumchi2=0;
3273       Float_t sum=0;
3274       for (Int_t i=0;i<6;i++){
3275         //
3276         Float_t c1 =0.; // conflict coeficients
3277         Float_t c2 =0.;
3278         if (clist1[i]&&clist2[i]){
3279           Float_t deltan = 0;
3280           if (i<2 || i>3){
3281             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
3282           }
3283           else{
3284             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
3285           }
3286           c1  = 2./TMath::Max(3.+deltan,2.);
3287           c2  = 2./TMath::Max(3.+deltan,2.);
3288         }
3289         else{
3290           if (clist1[i]){
3291             Float_t deltan = 0;
3292             if (i<2 || i>3){
3293               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
3294             }
3295             else{
3296               deltan = (clist1[i]->GetNz()-nz1[i]);
3297             }
3298             c1  = 2./TMath::Max(3.+deltan,2.);
3299             c2  = 0;
3300           }
3301
3302           if (clist2[i]){
3303             Float_t deltan = 0;
3304             if (i<2 || i>3){
3305               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
3306             }
3307             else{
3308               deltan = (clist2[i]->GetNz()-nz2[i]);
3309             }
3310             c2  = 2./TMath::Max(3.+deltan,2.);
3311             c1  = 0;
3312           }
3313         }
3314         //
3315         chi21=0;chi22=0;
3316         if (TMath::Abs(track1->GetDy(i))>0.) {
3317           chi21 = (track1->GetDy(i)/track1->GetSigmaY(i))*(track1->GetDy(i)/track1->GetSigmaY(i))+
3318             (track1->GetDz(i)/track1->GetSigmaZ(i))*(track1->GetDz(i)/track1->GetSigmaZ(i));
3319           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
3320           //  (track1->GetDz(i)*track1->GetDz(i))/(nerrz[i]*nerrz[i]);
3321         }else{
3322           if (TMath::Abs(track1->GetSigmaY(i)>0.)) c1=1;
3323         }
3324         //
3325         if (TMath::Abs(track2->GetDy(i))>0.) {
3326           chi22 = (track2->GetDy(i)/track2->GetSigmaY(i))*(track2->GetDy(i)/track2->GetSigmaY(i))+
3327             (track2->GetDz(i)/track2->GetSigmaZ(i))*(track2->GetDz(i)/track2->GetSigmaZ(i));
3328           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
3329           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
3330         }
3331         else{
3332           if (TMath::Abs(track2->GetSigmaY(i)>0.)) c2=1;
3333         }
3334         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
3335         if (chi21>0) sum+=w1;
3336         if (chi22>0) sum+=w2;
3337         conflict+=(c1+c2);
3338       }
3339       Double_t norm = sum-w1*track1->GetNSkipped()-w2*track2->GetNSkipped();
3340       if (norm<0) norm =1/(w1*track1->GetNSkipped()+w2*track2->GetNSkipped());
3341       Double_t normchi2 = 2*conflict+sumchi2/sum;
3342       if ( normchi2 <maxchi2 ){
3343         index1 = itrack1;
3344         index2 = itrack2;
3345         maxconflicts = conflict;
3346         maxchi2 = normchi2;
3347       }
3348     }
3349     UnRegisterClusterTracks(track1,trackID1);
3350   }
3351   //
3352   //  if (maxconflicts<4 && maxchi2<th0){
3353   if (maxchi2<th0*2.){
3354     Float_t orig = track10->GetFakeRatio()*track10->GetNumberOfClusters();
3355     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
3356     track1->SetChi2MIP(5,maxconflicts);
3357     track1->SetChi2MIP(6,maxchi2);
3358     track1->SetChi2MIP(7,0.01+orig-(track1->GetFakeRatio()*track1->GetNumberOfClusters()));
3359     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
3360     track1->SetChi2MIP(8,index1);
3361     fBestTrackIndex[trackID1] =index1;
3362     UpdateESDtrack(track1, AliESDtrack::kITSin);
3363     original->SetWinner(track1);
3364   }
3365   else if (track10->GetChi2MIP(0)<th1){
3366     track10->SetChi2MIP(5,maxconflicts);
3367     track10->SetChi2MIP(6,maxchi2);
3368     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
3369     UpdateESDtrack(track10,AliESDtrack::kITSin);
3370     original->SetWinner(track10);
3371   }
3372
3373   for (Int_t itrack=0;itrack<entries1;itrack++){
3374     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
3375     UnRegisterClusterTracks(track,trackID1);
3376   }
3377   //
3378   for (Int_t itrack=0;itrack<entries2;itrack++){
3379     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
3380     UnRegisterClusterTracks(track,trackID2);
3381   }
3382
3383   if (track10->GetConstrain()&&track10->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3384       &&track10->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3385     //  if (track10->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track10->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3386   //    &&track10->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track10->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3387     RegisterClusterTracks(track10,trackID1);
3388   }
3389   if (track20->GetConstrain()&&track20->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3390       &&track20->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3391     //if (track20->fChi2MIP[0]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track20->fChi2MIP[1]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3392     //  &&track20->fChi2MIP[2]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&track20->fChi2MIP[3]<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3393     RegisterClusterTracks(track20,trackID2);
3394   }
3395   return track10;
3396
3397 }
3398 //------------------------------------------------------------------------
3399 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
3400   //--------------------------------------------------------------------
3401   // This function marks clusters assigned to the track
3402   //--------------------------------------------------------------------
3403   AliTracker::UseClusters(t,from);
3404
3405   AliITSRecPoint *c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(0));
3406   //if (c->GetQ()>2) c->Use();
3407   if (c->GetSigmaZ2()>0.1) c->Use();
3408   c=(AliITSRecPoint *)GetCluster(t->GetClusterIndex(1));
3409   //if (c->GetQ()>2) c->Use();
3410   if (c->GetSigmaZ2()>0.1) c->Use();
3411
3412 }
3413 //------------------------------------------------------------------------
3414 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
3415 {
3416   //------------------------------------------------------------------
3417   // add track to the list of hypothesys
3418   //------------------------------------------------------------------
3419
3420   //
3421   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3422   if (!array) {
3423     array = new TObjArray(10);
3424     fTrackHypothesys.AddAtAndExpand(array,esdindex);
3425   }
3426   array->AddLast(track);
3427 }
3428 //------------------------------------------------------------------------
3429 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
3430 {
3431   //-------------------------------------------------------------------
3432   // compress array of track hypothesys
3433   // keep only maxsize best hypothesys
3434   //-------------------------------------------------------------------
3435   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
3436   if (! (fTrackHypothesys.At(esdindex)) ) return;
3437   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3438   Int_t entries = array->GetEntriesFast();
3439   //
3440   //- find preliminary besttrack as a reference
3441   Float_t minchi2=10000;
3442   Int_t maxn=0;
3443   AliITStrackMI * besttrack=0;
3444   //
3445   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
3446     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3447     if (!track) continue;
3448     Float_t chi2 = NormalizedChi2(track,0);
3449     //
3450     Int_t tpcLabel=track->GetESDtrack()->GetTPCLabel();
3451     track->SetLabel(tpcLabel);
3452     CookdEdx(track);
3453     track->SetFakeRatio(1.);
3454     CookLabel(track,0.); //For comparison only
3455     //
3456     //if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&track->fFakeRatio==0){
3457     if (chi2<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3458       if (track->GetNumberOfClusters()<maxn) continue;
3459       maxn = track->GetNumberOfClusters();
3460       //      if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2 *= track->GetChi2MIP(3); // RS
3461       if (chi2<minchi2){
3462         minchi2=chi2;
3463         besttrack=track;
3464       }
3465     }
3466     else{
3467       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3468         delete array->RemoveAt(itrack);
3469       }
3470     }
3471   }
3472   if (!besttrack) return;
3473   //
3474   //
3475   //take errors of best track as a reference
3476   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
3477   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3478   for (Int_t j=0;j<6;j++) {
3479     if (besttrack->GetClIndex(j)>=0){
3480       erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3481       errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3482       ny[j]   = besttrack->GetNy(j);
3483       nz[j]   = besttrack->GetNz(j);
3484     }
3485   }
3486   //
3487   // calculate normalized chi2
3488   //
3489   Float_t * chi2        = new Float_t[entries];
3490   Int_t * index         = new Int_t[entries];
3491   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3492   for (Int_t itrack=0;itrack<entries;itrack++){
3493     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3494     if (track){
3495       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3496       double chi2t = GetNormalizedChi2(track, mode);
3497       track->SetChi2MIP(0,chi2t);
3498       if (chi2t<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3499         if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3500         chi2[itrack] = chi2t;
3501       }
3502       else{
3503         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3504           delete array->RemoveAt(itrack);
3505         }
3506       }
3507     }
3508   }
3509   //
3510   TMath::Sort(entries,chi2,index,kFALSE);
3511   besttrack = (AliITStrackMI*)array->At(index[0]);
3512   if(besttrack) AliDebug(2,Form("ncls best track %d\n",besttrack->GetNumberOfClusters()));
3513   if (besttrack&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)){
3514     for (Int_t j=0;j<6;j++){
3515       if (besttrack->GetClIndex(j)>=0){
3516         erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3517         errz[j] = besttrack->GetSigmaZ(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3518         ny[j]   = besttrack->GetNy(j);
3519         nz[j]   = besttrack->GetNz(j);
3520       }
3521     }
3522   }
3523   //
3524   // calculate one more time with updated normalized errors
3525   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
3526   for (Int_t itrack=0;itrack<entries;itrack++){
3527     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3528     if (track){
3529       double chi2t = GetNormalizedChi2(track, mode);
3530       track->SetChi2MIP(0,chi2t);
3531       AliDebug(2,Form("track %d  ncls %d\n",itrack,track->GetNumberOfClusters()));
3532       if (track->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)) {
3533         if (fSelectBestMIP03 && track->GetChi2MIP(3)>0) chi2t *= track->GetChi2MIP(3); // RS
3534         chi2[itrack] = chi2t;  //-0*(track->GetNumberOfClusters()+track->GetNDeadZone());
3535       }
3536       else {
3537         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3538           delete array->RemoveAt(itrack);
3539         }
3540       }
3541     }
3542   }
3543   entries = array->GetEntriesFast();
3544   //
3545   //
3546   if (entries>0){
3547     TObjArray * newarray = new TObjArray();
3548     TMath::Sort(entries,chi2,index,kFALSE);
3549     besttrack = (AliITStrackMI*)array->At(index[0]);
3550     if (besttrack){
3551       AliDebug(2,Form("ncls best track %d     %f   %f\n",besttrack->GetNumberOfClusters(),besttrack->GetChi2MIP(0),chi2[index[0]]));
3552       //
3553       for (Int_t j=0;j<6;j++){
3554         if (besttrack->GetNz(j)>0&&besttrack->GetNy(j)>0){
3555           erry[j] = besttrack->GetSigmaY(j); erry[j+6] = besttrack->GetSigmaY(j+6);
3556           errz[j] = besttrack->GetSigmaZ(j); errz[j+6] = besttrack->GetSigmaZ(j+6);
3557           ny[j]   = besttrack->GetNy(j);
3558           nz[j]   = besttrack->GetNz(j);
3559         }
3560       }
3561       besttrack->SetChi2MIP(0,GetNormalizedChi2(besttrack,mode));
3562       minchi2 = TMath::Min(besttrack->GetChi2MIP(0)+5.+besttrack->GetNUsed(), double(AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)));
3563       Float_t minn = besttrack->GetNumberOfClusters()-3;
3564       Int_t accepted=0;
3565       for (Int_t i=0;i<entries;i++){
3566         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
3567         if (!track) continue;
3568         if (accepted>maxcut) break;
3569         track->SetChi2MIP(0,GetNormalizedChi2(track,mode));
3570         if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3571           if (track->GetNumberOfClusters()<6 && (track->GetChi2MIP(0)+track->GetNUsed()>minchi2)){
3572             delete array->RemoveAt(index[i]);
3573             continue;
3574           }
3575         }
3576         Bool_t shortbest = !track->GetConstrain() && track->GetNumberOfClusters()<6;
3577         if ((track->GetChi2MIP(0)+track->GetNUsed()<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
3578           if (!shortbest) accepted++;
3579           //
3580           newarray->AddLast(array->RemoveAt(index[i]));
3581           for (Int_t j=0;j<6;j++){
3582             if (nz[j]==0){
3583               erry[j] = track->GetSigmaY(j); erry[j+6] = track->GetSigmaY(j+6);
3584               errz[j] = track->GetSigmaZ(j); errz[j]   = track->GetSigmaZ(j+6);
3585               ny[j]   = track->GetNy(j);
3586               nz[j]   = track->GetNz(j);
3587             }
3588           }
3589         }
3590         else{
3591           delete array->RemoveAt(index[i]);
3592         }
3593       }
3594       array->Delete();
3595       delete fTrackHypothesys.RemoveAt(esdindex);
3596       fTrackHypothesys.AddAt(newarray,esdindex);
3597     }
3598     else{
3599       array->Delete();
3600       delete fTrackHypothesys.RemoveAt(esdindex);
3601     }
3602   }
3603   delete [] chi2;
3604   delete [] index;
3605 }
3606 //------------------------------------------------------------------------
3607 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
3608 {
3609   //-------------------------------------------------------------
3610   // try to find best hypothesy
3611   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
3612   // RS: optionally changing this to product of Chi2MIP(0)*Chi2MIP(3) == (chi2*chi2_interpolated)
3613   //-------------------------------------------------------------
3614   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
3615   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
3616   if (!array) return 0;
3617   Int_t entries = array->GetEntriesFast();
3618   if (!entries) return 0;
3619   Float_t minchi2 = 100000;
3620   AliITStrackMI * besttrack=0;
3621   //
3622   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
3623   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
3624   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
3625   Double_t ersVtx[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
3626   //
3627   for (Int_t i=0;i<entries;i++){
3628     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3629     if (!track) continue;
3630     Float_t sigmarfi,sigmaz;
3631     GetDCASigma(track,sigmarfi,sigmaz);
3632     track->SetDnorm(0,sigmarfi);
3633     track->SetDnorm(1,sigmaz);
3634     //
3635     track->SetChi2MIP(1,1000000);
3636     track->SetChi2MIP(2,1000000);
3637     track->SetChi2MIP(3,1000000);
3638     //
3639     // backtrack
3640     backtrack = new(backtrack) AliITStrackMI(*track);
3641     if (track->GetConstrain()) {
3642       if (!CorrectForPipeMaterial(backtrack,"inward")) continue;
3643       if (AliITSReconstructor::GetRecoParam()->GetImproveWithVertex()) {
3644         if (fUseImproveKalman) {if (!backtrack->ImproveKalman(xyzVtx,ersVtx,0,0,0)) continue;}
3645         else                   {if (!backtrack->Improve(0,xyzVtx,ersVtx)) continue;}
3646       }
3647       backtrack->ResetCovariance(10.);
3648     }else{
3649       backtrack->ResetCovariance(10.);
3650     }
3651     backtrack->ResetClusters();
3652
3653     Double_t x = original->GetX();
3654     if (!RefitAt(x,backtrack,track)) continue;
3655     //
3656     track->SetChi2MIP(1,NormalizedChi2(backtrack,0));
3657     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
3658     if (track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)*6.)  continue;
3659     track->SetChi22(GetMatchingChi2(backtrack,original));
3660
3661     if ((track->GetConstrain()) && track->GetChi22()>90.)  continue;
3662     if ((!track->GetConstrain()) && track->GetChi22()>30.)  continue;
3663     if ( track->GetChi22()/track->GetNumberOfClusters()>11.)  continue;
3664
3665
3666     if  (!(track->GetConstrain())&&track->GetChi2MIP(1)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1))  continue;
3667     //
3668     //forward track - without constraint
3669     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
3670     forwardtrack->ResetClusters();
3671     x = track->GetX();
3672     if (!RefitAt(x,forwardtrack,track) && fSelectBestMIP03) continue;  // w/o fwd track MIP03 is meaningless
3673     track->SetChi2MIP(2,NormalizedChi2(forwardtrack,0));
3674     if  (track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)*6.0)  continue;
3675     if  (!(track->GetConstrain())&&track->GetChi2MIP(2)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2))  continue;
3676
3677     //track->fD[0] = forwardtrack->GetD(GetX(),GetY());
3678     //track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3679     forwardtrack->GetDZ(GetX(),GetY(),GetZ(),track->GetDP());   //I.B.
3680     forwardtrack->SetD(0,track->GetD(0));
3681     forwardtrack->SetD(1,track->GetD(1));
3682     {
3683       Int_t list[6];
3684       AliITSRecPoint* clist[6];
3685       track->SetChi2MIP(4,GetNumberOfSharedClusters(track,esdindex,list,clist));
3686       if ( (!track->GetConstrain()) && track->GetChi2MIP(4)>1.0) continue;
3687     }
3688
3689     track->SetChi2MIP(3,GetInterpolatedChi2(forwardtrack,backtrack));
3690     if  ( (track->GetChi2MIP(3)>6.*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) continue;
3691     if  ( (!track->GetConstrain()) && (track->GetChi2MIP(3)>2*AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3))) {
3692       track->SetChi2MIP(3,1000);
3693       continue;
3694     }
3695     Double_t chi2 = track->GetChi2MIP(0); // +track->GetNUsed();    //RS
3696     if (fSelectBestMIP03) chi2 *= track->GetChi2MIP(3);
3697     else chi2 += track->GetNUsed();
3698     //
3699     for (Int_t ichi=0;ichi<5;ichi++){
3700       forwardtrack->SetChi2MIP(ichi, track->GetChi2MIP(ichi));
3701     }
3702     if (chi2 < minchi2){
3703       //besttrack = new AliITStrackMI(*forwardtrack);
3704       besttrack = track;
3705       besttrack->SetLabel(track->GetLabel());
3706       besttrack->SetFakeRatio(track->GetFakeRatio());
3707       minchi2   = chi2;
3708       //original->fD[0] = forwardtrack->GetD(GetX(),GetY());
3709       //original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
3710       forwardtrack->GetDZ(GetX(),GetY(),GetZ(),original->GetDP());    //I.B.
3711     }
3712   }
3713   delete backtrack;
3714   delete forwardtrack;
3715
3716   if (!besttrack)  return 0;
3717
3718   Int_t accepted=0;
3719   for (Int_t i=0;i<entries;i++){
3720     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3721
3722     if (!track) continue;
3723
3724     if (accepted>checkmax || track->GetChi2MIP(3)>AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)*6. ||
3725         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)
3726         // RS: don't apply this cut when fSelectBestMIP03 is on
3727         || (!fSelectBestMIP03 && (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*besttrack->GetNUsed()+3.))
3728         ){
3729       if (track->GetConstrain() || track->GetNumberOfClusters()>5){  //keep best short tracks - without vertex constrain
3730         delete array->RemoveAt(i);
3731         continue;
3732       }
3733     }
3734     else{
3735       accepted++;
3736     }
3737   }
3738   //
3739   array->Compress();
3740   SortTrackHypothesys(esdindex,checkmax,1);
3741
3742   array = (TObjArray*) fTrackHypothesys.At(esdindex);
3743   if (!array) return 0; // PH What can be the reason? Check SortTrackHypothesys
3744   besttrack = (AliITStrackMI*)array->At(0);
3745   if (!besttrack)  return 0;
3746   besttrack->SetChi2MIP(8,0);
3747   fBestTrackIndex[esdindex]=0;
3748   entries = array->GetEntriesFast();
3749   AliITStrackMI *longtrack =0;
3750   minchi2 =1000;
3751   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->GetNDeadZone();
3752   for (Int_t itrack=entries-1;itrack>0;itrack--) {
3753     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
3754     if (!track->GetConstrain()) continue;
3755     if (track->GetNumberOfClusters()+track->GetNDeadZone()<minn) continue;
3756     if (track->GetChi2MIP(0)-besttrack->GetChi2MIP(0)>0.0) continue;
3757     if (track->GetChi2MIP(0)>4.) continue;
3758     minn = track->GetNumberOfClusters()+track->GetNDeadZone();
3759     longtrack =track;
3760   }
3761   //if (longtrack) besttrack=longtrack;
3762   //
3763   // RS do shared cluster analysis here only if the new sharing analysis is not requested
3764   //RRR if (fFlagFakes) return besttrack;
3765
3766   Int_t list[6];
3767   AliITSRecPoint * clist[6];
3768   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3769   if (besttrack->GetConstrain()&&besttrack->GetChi2MIP(0)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(0)&&besttrack->GetChi2MIP(1)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(1)
3770       &&besttrack->GetChi2MIP(2)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(2)&&besttrack->GetChi2MIP(3)<AliITSReconstructor::GetRecoParam()->GetMaxChi2PerCluster(3)){
3771     RegisterClusterTracks(besttrack,esdindex);
3772   }
3773   //
3774   //
3775   if (shared>0.0){
3776     Int_t nshared;
3777     Int_t overlist[6];
3778     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
3779     if (sharedtrack>=0){
3780       //
3781       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5,original);
3782       if (besttrack){
3783         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
3784       }
3785       else return 0;
3786     }
3787   }
3788
3789   if (shared>2.5) return 0;
3790   if (shared>1.0) return besttrack;
3791   //
3792   // Don't sign clusters if not gold track
3793   //
3794   if (!besttrack->IsGoldPrimary()) return besttrack;
3795   if (besttrack->GetESDtrack()->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
3796   //
3797   if (fConstraint[fPass]){
3798     //
3799     // sign clusters
3800     //
3801     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
3802     for (Int_t i=0;i<6;i++){
3803       Int_t index = besttrack->GetClIndex(i);
3804       if (index<0) continue;
3805       Int_t ilayer =  (index & 0xf0000000) >> 28;
3806       if (besttrack->GetSigmaY(ilayer)<0.00000000001) continue;
3807       AliITSRecPoint *c = (AliITSRecPoint*)GetCluster(index);
3808       if (!c) continue;
3809       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
3810       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
3811       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
3812       if ( ilayer>2&& AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(ilayer))
3813         if (besttrack->GetNormQ(ilayer)/besttrack->GetExpQ()>1.5) continue;
3814       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
3815
3816       Bool_t cansign = kTRUE;
3817       for (Int_t itrack=0;itrack<entries; itrack++){
3818         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
3819         if (!track) continue;
3820         if (track->GetChi2MIP(0)>besttrack->GetChi2MIP(0)+2.*shared+1.) break;
3821         if ( (track->GetClIndex(ilayer)>=0) && (track->GetClIndex(ilayer)!=besttrack->GetClIndex(ilayer))){
3822           cansign = kFALSE;
3823           break;
3824         }
3825       }
3826       if (cansign){
3827         if (TMath::Abs(besttrack->GetDy(ilayer)/besttrack->GetSigmaY(ilayer))>3.) continue;
3828         if (TMath::Abs(besttrack->GetDz(ilayer)/besttrack->GetSigmaZ(ilayer))>3.) continue;
3829         if (!c->IsUsed()) c->Use();
3830       }
3831     }
3832   }
3833   return besttrack;
3834 }
3835 //------------------------------------------------------------------------
3836 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
3837 {
3838   //
3839   // get "best" hypothesys
3840   //
3841
3842   Int_t nentries = itsTracks.GetEntriesFast();
3843   for (Int_t i=0;i<nentries;i++){
3844     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
3845     if (!track) continue;
3846     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
3847     if (!array) continue;
3848     if (array->GetEntriesFast()<=0) continue;
3849     //
3850     AliITStrackMI* longtrack=0;
3851     Float_t minn=0;
3852     Float_t maxchi2=1000;
3853     for (Int_t j=0;j<array->GetEntriesFast();j++){
3854       AliITStrackMI* trackHyp = (AliITStrackMI*)array->At(j);
3855       if (!trackHyp) continue;
3856       if (trackHyp->GetGoldV0()) {
3857         longtrack = trackHyp;   //gold V0 track taken
3858         break;
3859       }
3860       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()<minn) continue;
3861       Float_t chi2 = trackHyp->GetChi2MIP(0);
3862       if (fSelectBestMIP03) chi2 *= trackHyp->GetChi2MIP(3);
3863       if (trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone()>minn) maxchi2 = chi2; //trackHyp->GetChi2MIP(0);
3864       //
3865       if (fAfterV0){ // ??? RS
3866         if (!trackHyp->GetGoldV0()&&trackHyp->GetConstrain()==kFALSE) chi2+=5;
3867       }
3868       if (chi2 > maxchi2) continue;
3869       minn = trackHyp->GetNumberOfClusters()+trackHyp->GetNDeadZone();
3870       if (fSelectBestMIP03) minn++; // allow next to longest to win
3871       maxchi2 = chi2;
3872       longtrack=trackHyp;
3873     }
3874     //
3875     //
3876     //
3877     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
3878     if (!longtrack) {longtrack = besttrack;}
3879     else besttrack= longtrack;
3880     //
3881     if (besttrack) {
3882       Int_t list[6];
3883       AliITSRecPoint * clist[6];
3884       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
3885       //
3886       track->SetNUsed(shared);
3887       track->SetNSkipped(besttrack->GetNSkipped());
3888       track->SetChi2MIP(0,besttrack->GetChi2MIP(0));
3889       if (shared>0) {
3890         if(!AliITSReconstructor::GetRecoParam()->GetAllowSharedClusters()) continue;
3891         Int_t nshared;
3892         Int_t overlist[6];
3893         //
3894         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
3895         //if (sharedtrack==-1) sharedtrack=0;
3896         if (sharedtrack>=0) {
3897           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5,track);
3898         }
3899       }
3900       if (besttrack&&fAfterV0) {
3901         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3902         track->SetWinner(besttrack);
3903       }
3904       if (besttrack) {
3905         if (fConstraint[fPass]) {
3906           UpdateESDtrack(besttrack,AliESDtrack::kITSin);
3907           track->SetWinner(besttrack);
3908         }
3909         if (besttrack->GetChi2MIP(0)+besttrack->GetNUsed()>1.5 && fConstraint[fPass]) {
3910           if ( TMath::Abs(besttrack->GetD(0))>0.1 ||
3911                TMath::Abs(besttrack->GetD(1))>0.1 ) track->SetReconstructed(kFALSE);
3912         }
3913       }
3914     }
3915   }
3916 }
3917
3918 //------------------------------------------------------------------------
3919 void AliITStrackerMI::FlagFakes(const TObjArray &itsTracks)
3920 {
3921   //
3922   // RS: flag those tracks which are suxpected to have fake clusters
3923   //
3924   const double kThreshPt = 0.5;
3925   AliRefArray *refArr[6];
3926   //
3927   for (int i=0;i<6;i++) {
3928     int ncl = fgLayers[i].GetNumberOfClusters();
3929     refArr[i] = new AliRefArray(ncl,TMath::Min(ncl,1000));
3930   }
3931   Int_t nentries = itsTracks.GetEntriesFast();
3932   //
3933   // fill cluster->track associations
3934   for (Int_t itr=0;itr<nentries;itr++){
3935     AliITStrackMI* track = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3936     if (!track) continue;
3937     AliITStrackMI* trackITS = track->GetWinner();
3938     if (!trackITS) continue;
3939     for (int il=trackITS->GetNumberOfClusters();il--;) {
3940       int idx = trackITS->GetClusterIndex(il);
3941       Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3942       //      if (c>fgLayers[l].GetNumberOfClusters()) continue;
3943       refArr[l]->AddReference(c, itr);
3944     }
3945   }
3946   //
3947   const UInt_t kMaxRef = 100;
3948   UInt_t crefs[kMaxRef];
3949   Int_t ncrefs=0;
3950   // process tracks with shared clusters
3951   for (int itr=0;itr<nentries;itr++){
3952     AliITStrackMI* track0 = (AliITStrackMI*)itsTracks.UncheckedAt(itr);
3953     AliITStrackMI* trackH0 = track0->GetWinner();
3954     if (!trackH0) continue;
3955     AliESDtrack* esd0 = track0->GetESDtrack();
3956     //
3957     for (int il=0;il<trackH0->GetNumberOfClusters();il++) {
3958       int idx = trackH0->GetClusterIndex(il);
3959       Int_t l=(idx & 0xf0000000) >> 28, c=(idx & 0x0fffffff) >> 00;
3960       ncrefs = refArr[l]->GetReferences(c,crefs,kMaxRef);
3961       if (ncrefs<2) continue; // there will be always self-reference, for sharing needs at least 2
3962       esd0->SetITSSharedFlag(l);
3963       for (int ir=ncrefs;ir--;) {
3964         if (int(crefs[ir]) <= itr) continue; // ==:selfreference, <: the same pair will be checked with >
3965         AliITStrackMI* track1 = (AliITStrackMI*)itsTracks.UncheckedAt(crefs[ir]);
3966         AliITStrackMI* trackH1 = track1->GetWinner();
3967         AliESDtrack* esd1 = track1->GetESDtrack();
3968         esd1->SetITSSharedFlag(l);
3969         //
3970         double pt0 = trackH0->Pt(), pt1 = trackH1->Pt(), res = 0.;
3971         if      (pt0>kThreshPt && pt0-pt1>0.2+0.2*(pt0-kThreshPt) ) res = -100;
3972         else if (pt1>kThreshPt && pt1-pt0>0.2+0.2*(pt1-kThreshPt) ) res = 100;
3973
3974         // select the one with smallest chi2's product
3975         res += trackH0->GetChi2MIP(0)*trackH0->GetChi2MIP(3);
3976         res -= trackH1->GetChi2MIP(0)*trackH1->GetChi2MIP(3);
3977         //
3978         if (res<0) esd1->SetITSFakeFlag();  // esd0 is winner
3979         else       esd0->SetITSFakeFlag();  // esd1 is winner
3980       }
3981       //
3982     }
3983     //
3984   }
3985   //
3986   for (int i=6;i--;) delete refArr[i];
3987 }
3988
3989 //------------------------------------------------------------------------
3990 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
3991   //--------------------------------------------------------------------
3992   //This function "cooks" a track label. If label<0, this track is fake.
3993   //--------------------------------------------------------------------
3994   Int_t tpcLabel=-1;
3995
3996   if (track->GetESDtrack()){
3997     tpcLabel = track->GetESDtrack()->GetTPCLabel();
3998     ULong_t trStatus=track->GetESDtrack()->GetStatus();
3999     if(!(trStatus&AliESDtrack::kTPCin)) tpcLabel=track->GetLabel(); // for ITSsa tracks
4000   }
4001    track->SetChi2MIP(9,0);
4002    Int_t nwrong=0;
4003    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4004      Int_t cindex = track->GetClusterIndex(i);
4005      Int_t l=(cindex & 0xf0000000) >> 28;
4006      AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4007      Int_t isWrong=1;
4008      for (Int_t ind=0;ind<3;ind++){
4009        if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
4010        //AliDebug(2,Form("icl %d  ilab %d lab %d",i,ind,cl->GetLabel(ind)));
4011      }
4012      track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4013      nwrong+=isWrong;
4014    }
4015    Int_t nclusters = track->GetNumberOfClusters();
4016    if (nclusters > 0) //PH Some tracks don't have any cluster
4017      track->SetFakeRatio(double(nwrong)/double(nclusters));
4018    if (tpcLabel>0 && track->GetFakeRatio()>wrong) {
4019      track->SetLabel(-tpcLabel);
4020    } else {
4021      track->SetLabel(tpcLabel);
4022    }
4023    AliDebug(2,Form(" nls %d wrong %d  label %d  tpcLabel %d\n",nclusters,nwrong,track->GetLabel(),tpcLabel));
4024
4025 }
4026 //------------------------------------------------------------------------
4027 void AliITStrackerMI::CookdEdx(AliITStrackMI* track){
4028   //
4029   // Fill the dE/dx in this track
4030   //
4031   track->SetChi2MIP(9,0);
4032   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
4033     Int_t cindex = track->GetClusterIndex(i);
4034     Int_t l=(cindex & 0xf0000000) >> 28;
4035     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
4036     Int_t lab = TMath::Abs(track->GetESDtrack()->GetTPCLabel());
4037     Int_t isWrong=1;
4038     for (Int_t ind=0;ind<3;ind++){
4039       if (cl->GetLabel(ind)==lab) isWrong=0;
4040     }
4041     track->SetChi2MIP(9,track->GetChi2MIP(9)+isWrong*(2<<l));
4042   }
4043   Double_t low=0.;
4044   Double_t up=0.51;
4045   track->CookdEdx(low,up);
4046 }
4047 //------------------------------------------------------------------------
4048 void AliITStrackerMI::MakeCoefficients(Int_t ntracks){
4049   //
4050   // Create some arrays
4051   //
4052   if (fCoefficients) delete []fCoefficients;
4053   fCoefficients = new Float_t[ntracks*48];
4054   for (Int_t i=0;i<ntracks*48;i++) fCoefficients[i]=-1.;
4055 }
4056 //------------------------------------------------------------------------
4057 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4058 {
4059   //
4060   // Compute predicted chi2
4061   //
4062   // Take into account the mis-alignment (bring track to cluster plane)
4063   Double_t xTrOrig=track->GetX();
4064   if (!track->Propagate(xTrOrig+cluster->GetX())) return 1000.;
4065   Float_t erry,errz,covyz;
4066   Float_t theta = track->GetTgl();
4067   Float_t phi   = track->GetSnp();
4068   phi *= TMath::Sqrt(1./((1.-phi)*(1.+phi)));
4069   AliITSClusterParam::GetError(layer,cluster,theta,phi,track->GetExpQ(),erry,errz,covyz);
4070   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
4071   AliDebug(2,Form(" chi2: tr-cl   %f  %f   tr X %f cl X %f",track->GetY()-cluster->GetY(),track->GetZ()-cluster->GetZ(),track->GetX(),cluster->GetX()));
4072   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz,covyz);
4073   // Bring the track back to detector plane in ideal geometry
4074   // [mis-alignment will be accounted for in UpdateMI()]
4075   if (!track->Propagate(xTrOrig)) return 1000.;
4076   Float_t ny,nz;
4077   AliITSClusterParam::GetNTeor(layer,cluster,theta,phi,ny,nz);
4078   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
4079   if (delta>1){
4080     chi2+=0.5*TMath::Min(delta/2,2.);
4081     chi2+=2.*cluster->GetDeltaProbability();
4082   }
4083   //
4084   track->SetNy(layer,ny);
4085   track->SetNz(layer,nz);
4086   track->SetSigmaY(layer,erry);
4087   track->SetSigmaZ(layer, errz);
4088   track->SetSigmaYZ(layer,covyz);
4089   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
4090   track->SetNormQ(layer,cluster->GetQ()/TMath::Sqrt((1.+ track->GetTgl()*track->GetTgl())/((1.-track->GetSnp())*(1.+track->GetSnp()))));
4091   return chi2;
4092
4093 }
4094 //------------------------------------------------------------------------
4095 Int_t AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSRecPoint* cl,Double_t chi2,Int_t index) const
4096 {
4097   //
4098   // Update ITS track
4099   //
4100   Int_t layer = (index & 0xf0000000) >> 28;
4101   track->SetClIndex(layer, index);
4102   if (layer>1&&AliITSReconstructor::GetRecoParam()->GetUseAmplitudeInfo(layer)) {
4103     if (track->GetNormQ(layer)/track->GetExpQ()<0.5 ) {
4104       chi2+= (0.5-track->GetNormQ(layer)/track->GetExpQ())*10.;
4105       track->SetdEdxMismatch(track->GetdEdxMismatch()+(0.5-track->GetNormQ(layer)/track->GetExpQ())*10.);
4106     }
4107   }
4108
4109   if (TMath::Abs(cl->GetQ())<1.e-13) return 0;  // ingore the "virtual" clusters
4110
4111
4112   // Take into account the mis-alignment (bring track to cluster plane)
4113   Double_t xTrOrig=track->GetX();
4114   Float_t clxyz[3]; cl->GetGlobalXYZ(clxyz);Double_t trxyz[3]; track->GetXYZ(trxyz);
4115   AliDebug(2,Form("gtr %f %f %f",trxyz[0],trxyz[1],trxyz[2]));
4116   AliDebug(2,Form("gcl %f %f %f",clxyz[0],clxyz[1],clxyz[2]));
4117   AliDebug(2,Form(" xtr %f  xcl %f",track->GetX(),cl->GetX()));
4118
4119   if (!track->Propagate(xTrOrig+cl->GetX())) return 0;
4120
4121   AliCluster c(*cl);
4122   c.SetSigmaY2(track->GetSigmaY(layer)*track->GetSigmaY(layer));
4123   c.SetSigmaZ2(track->GetSigmaZ(layer)*track->GetSigmaZ(layer));
4124   c.SetSigmaYZ(track->GetSigmaYZ(layer));
4125
4126
4127   Int_t updated = track->UpdateMI(&c,chi2,index);
4128   // Bring the track back to detector plane in ideal geometry
4129   if (!track->Propagate(xTrOrig)) return 0;
4130
4131   if(!updated) AliDebug(2,"update failed");
4132   return updated;
4133 }
4134
4135 //------------------------------------------------------------------------
4136 void AliITStrackerMI::GetDCASigma(const AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
4137 {
4138   //
4139   //DCA sigmas parameterization
4140   //to be paramterized using external parameters in future
4141   //
4142   //
4143   Double_t curv=track->GetC();
4144   sigmarfi = 0.0040+1.4 *TMath::Abs(curv)+332.*curv*curv;
4145   sigmaz   = 0.0110+4.37*TMath::Abs(curv);
4146 }
4147 //------------------------------------------------------------------------
4148 void AliITStrackerMI::SignDeltas(const TObjArray *clusterArray, Float_t vz)
4149 {
4150   //
4151   // Clusters from delta electrons?
4152   //
4153   Int_t entries = clusterArray->GetEntriesFast();
4154   if (entries<4) return;
4155   AliITSRecPoint* cluster = (AliITSRecPoint*)clusterArray->At(0);
4156   Int_t layer = cluster->GetLayer();
4157   if (layer>1) return;
4158   Int_t index[10000];
4159   Int_t ncandidates=0;
4160   Float_t r = (layer>0)? 7:4;
4161   //
4162   for (Int_t i=0;i<entries;i++){
4163     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(i);
4164     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
4165     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
4166     index[ncandidates] = i;  //candidate to belong to delta electron track
4167     ncandidates++;
4168     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
4169       cl0->SetDeltaProbability(1);
4170     }
4171   }
4172   //
4173   //
4174   //
4175   for (Int_t i=0;i<ncandidates;i++){
4176     AliITSRecPoint* cl0 = (AliITSRecPoint*)clusterArray->At(index[i]);
4177     if (cl0->GetDeltaProbability()>0.8) continue;
4178     //
4179     Int_t ncl = 0;
4180     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
4181     sumy=sumz=sumy2=sumyz=sumw=0.0;
4182     for (Int_t j=0;j<ncandidates;j++){
4183       if (i==j) continue;
4184       AliITSRecPoint* cl1 = (AliITSRecPoint*)clusterArray->At(index[j]);
4185       //
4186       Float_t dz = cl0->GetZ()-cl1->GetZ();
4187       Float_t dy = cl0->GetY()-cl1->GetY();
4188       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
4189         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
4190         y[ncl] = cl1->GetY();
4191         z[ncl] = cl1->GetZ();
4192         sumy+= y[ncl]*weight;
4193         sumz+= z[ncl]*weight;
4194         sumy2+=y[ncl]*y[ncl]*weight;
4195         sumyz+=y[ncl]*z[ncl]*weight;
4196         sumw+=weight;
4197         ncl++;
4198       }
4199     }
4200     if (ncl<4) continue;
4201     Float_t det = sumw*sumy2  - sumy*sumy;
4202     Float_t delta=1000;
4203     if (TMath::Abs(det)>0.01){
4204       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
4205       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
4206       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
4207     }
4208     else{
4209       Float_t z0  = sumyz/sumy;
4210       delta = TMath::Abs(cl0->GetZ()-z0);
4211     }
4212     if ( delta<0.05) {
4213       cl0->SetDeltaProbability(1-20.*delta);
4214     }
4215   }
4216 }
4217 //------------------------------------------------------------------------
4218 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
4219 {
4220   //
4221   // Update ESD track
4222   //
4223   track->UpdateESDtrack(flags);
4224   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->GetESDtrack()->GetITStrack());
4225   if (oldtrack) delete oldtrack;
4226   track->GetESDtrack()->SetITStrack(new AliITStrackMI(*track));
4227   // if (TMath::Abs(track->GetDnorm(1))<0.000000001){
4228   //   printf("Problem\n");
4229   // }
4230 }
4231 //------------------------------------------------------------------------
4232 Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
4233   //
4234   // Get nearest upper layer close to the point xr.
4235   // rough approximation
4236   //
4237   const Float_t kRadiuses[6]={4,6.5,15.03,24.,38.5,43.7};
4238   Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
4239   Int_t res =6;
4240   for (Int_t i=0;i<6;i++){
4241     if (radius<kRadiuses[i]){
4242       res =i;
4243       break;
4244     }
4245   }
4246   return res;
4247 }
4248 //------------------------------------------------------------------------
4249 void AliITStrackerMI::BuildMaterialLUT(TString material) {
4250   //--------------------------------------------------------------------
4251   // Fill a look-up table with mean material
4252   //--------------------------------------------------------------------
4253
4254   Int_t n=1000;
4255   Double_t mparam[7];
4256   Double_t point1[3],point2[3];
4257   Double_t phi,cosphi,sinphi,z;
4258   // 0-5 layers, 6 pipe, 7-8 shields
4259   Double_t rmin[9]={ 3.5, 5.5,13.0,22.0,35.0,41.0, 2.0, 8.0,25.0};
4260   Double_t rmax[9]={ 5.5, 8.0,17.0,26.0,41.0,47.0, 3.0,10.5,30.0};
4261
4262   Int_t ifirst=0,ilast=0;
4263   if(material.Contains("Pipe")) {
4264     ifirst=6; ilast=6;
4265   } else if(material.Contains("Shields")) {
4266     ifirst=7; ilast=8;
4267   } else if(material.Contains("Layers")) {
4268     ifirst=0; ilast=5;
4269   } else {
4270     Error("BuildMaterialLUT","Wrong layer name\n");
4271   }
4272
4273   for(Int_t imat=ifirst; imat<=ilast; imat++) {
4274     Double_t param[5]={0.,0.,0.,0.,0.};
4275     for (Int_t i=0; i<n; i++) {
4276       phi = 2.*TMath::Pi()*gRandom->Rndm();
4277       cosphi = TMath::Cos(phi); sinphi = TMath::Sin(phi);
4278       z = 14.*(-1.+2.*gRandom->Rndm()); // SPD barrel
4279       point1[0] = rmin[imat]*cosphi;
4280       point1[1] = rmin[imat]*sinphi;
4281       point1[2] = z;
4282       point2[0] = rmax[imat]*cosphi;
4283       point2[1] = rmax[imat]*sinphi;
4284       point2[2] = z;
4285       AliTracker::MeanMaterialBudget(point1,point2,mparam);
4286       for(Int_t j=0;j<5;j++) param[j]+=mparam[j];
4287     }
4288     for(Int_t j=0;j<5;j++) param[j]/=(Float_t)n;
4289     if(imat<=5) {
4290       fxOverX0Layer[imat] = param[1];
4291       fxTimesRhoLayer[imat] = param[0]*param[4];
4292     } else if(imat==6) {
4293       fxOverX0Pipe = param[1];
4294       fxTimesRhoPipe = param[0]*param[4];
4295     } else if(imat==7) {
4296       fxOverX0Shield[0] = param[1];
4297       fxTimesRhoShield[0] = param[0]*param[4];
4298     } else if(imat==8) {
4299       fxOverX0Shield[1] = param[1];
4300       fxTimesRhoShield[1] = param[0]*param[4];
4301     }
4302   }
4303   /*
4304   printf("%s\n",material.Data());
4305   printf("%f  %f\n",fxOverX0Pipe,fxTimesRhoPipe);
4306   printf("%f  %f\n",fxOverX0Shield[0],fxTimesRhoShield[0]);
4307   printf("%f  %f\n",fxOverX0Shield[1],fxTimesRhoShield[1]);
4308   printf("%f  %f\n",fxOverX0Layer[0],fxTimesRhoLayer[0]);
4309   printf("%f  %f\n",fxOverX0Layer[1],fxTimesRhoLayer[1]);
4310   printf("%f  %f\n",fxOverX0Layer[2],fxTimesRhoLayer[2]);
4311   printf("%f  %f\n",fxOverX0Layer[3],fxTimesRhoLayer[3]);
4312   printf("%f  %f\n",fxOverX0Layer[4],fxTimesRhoLayer[4]);
4313   printf("%f  %f\n",fxOverX0Layer[5],fxTimesRhoLayer[5]);
4314   */
4315   return;
4316 }
4317 //------------------------------------------------------------------------
4318 Int_t AliITStrackerMI::CorrectForPipeMaterial(AliITStrackMI *t,
4319                                               TString direction) {
4320   //-------------------------------------------------------------------
4321   // Propagate beyond beam pipe and correct for material
4322   // (material budget in different ways according to fUseTGeo value)
4323   // Add time if going outward (PropagateTo or PropagateToTGeo)
4324   //-------------------------------------------------------------------
4325
4326   // Define budget mode:
4327   // 0: material from AliITSRecoParam (hard coded)
4328   // 1: material from TGeo in one step (on the fly)
4329   // 2: material from lut
4330   // 3: material from TGeo in one step (same for all hypotheses)
4331   Int_t mode;
4332   switch(fUseTGeo) {
4333   case 0:
4334     mode=0;
4335     break;
4336   case 1:
4337     mode=1;
4338     break;
4339   case 2:
4340     mode=2;
4341     break;
4342   case 3:
4343     if(fTrackingPhase.Contains("Clusters2Tracks"))
4344       { mode=3; } else { mode=1; }
4345     break;
4346   case 4:
4347     if(fTrackingPhase.Contains("Clusters2Tracks"))
4348       { mode=3; } else { mode=2; }
4349     break;
4350   default:
4351     mode=0;
4352     break;
4353   }
4354   if(fTrackingPhase.Contains("Default")) mode=0;
4355
4356   Int_t index=fCurrentEsdTrack;
4357
4358   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4359   Double_t rToGo=(dir>0 ? AliITSRecoParam::GetrInsidePipe() : AliITSRecoParam::GetrOutsidePipe());
4360   Double_t xToGo;
4361   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4362
4363   Double_t xOverX0,x0,lengthTimesMeanDensity;
4364
4365   switch(mode) {
4366   case 0:
4367     xOverX0 = AliITSRecoParam::GetdPipe();
4368     x0 = AliITSRecoParam::GetX0Be();
4369     lengthTimesMeanDensity = xOverX0*x0;
4370     lengthTimesMeanDensity *= dir;
4371     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4372     break;
4373   case 1:
4374     if (!t->PropagateToTGeo(xToGo,1)) return 0;
4375     break;
4376   case 2:
4377     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
4378     xOverX0 = fxOverX0Pipe;
4379     lengthTimesMeanDensity = fxTimesRhoPipe;
4380     lengthTimesMeanDensity *= dir;
4381     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4382     break;
4383   case 3:
4384     if(!fxOverX0PipeTrks || index<0 || index>=fNtracks) Error("CorrectForPipeMaterial","Incorrect usage of UseTGeo option!\n");
4385     if(fxOverX0PipeTrks[index]<0) {
4386       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4387       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4388                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4389       fxOverX0PipeTrks[index] = TMath::Abs(xOverX0)/angle;
4390       fxTimesRhoPipeTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4391       return 1;
4392     }
4393     xOverX0 = fxOverX0PipeTrks[index];
4394     lengthTimesMeanDensity = fxTimesRhoPipeTrks[index];
4395     lengthTimesMeanDensity *= dir;
4396     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4397     break;
4398   }
4399
4400   return 1;
4401 }
4402 //------------------------------------------------------------------------
4403 Int_t AliITStrackerMI::CorrectForShieldMaterial(AliITStrackMI *t,
4404                                                 TString shield,
4405                                                 TString direction) {
4406   //-------------------------------------------------------------------
4407   // Propagate beyond SPD or SDD shield and correct for material
4408   // (material budget in different ways according to fUseTGeo value)
4409   // Add time if going outward (PropagateTo or PropagateToTGeo)
4410   //-------------------------------------------------------------------
4411
4412   // Define budget mode:
4413   // 0: material from AliITSRecoParam (hard coded)
4414   // 1: material from TGeo in steps of X cm (on the fly)
4415   //    X = AliITSRecoParam::GetStepSizeTGeo()
4416   // 2: material from lut
4417   // 3: material from TGeo in one step (same for all hypotheses)
4418   Int_t mode;
4419   switch(fUseTGeo) {
4420   case 0:
4421     mode=0;
4422     break;
4423   case 1:
4424     mode=1;
4425     break;
4426   case 2:
4427     mode=2;
4428     break;
4429   case 3:
4430     if(fTrackingPhase.Contains("Clusters2Tracks"))
4431       { mode=3; } else { mode=1; }
4432     break;
4433   case 4:
4434     if(fTrackingPhase.Contains("Clusters2Tracks"))
4435       { mode=3; } else { mode=2; }
4436     break;
4437   default:
4438     mode=0;
4439     break;
4440   }
4441   if(fTrackingPhase.Contains("Default")) mode=0;
4442
4443   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4444   Double_t rToGo;
4445   Int_t    shieldindex=0;
4446   if (shield.Contains("SDD")) { // SDDouter
4447     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(1) : AliITSRecoParam::GetrOutsideShield(1));
4448     shieldindex=1;
4449   } else if (shield.Contains("SPD")) {        // SPDouter
4450     rToGo=(dir>0 ? AliITSRecoParam::GetrInsideShield(0) : AliITSRecoParam::GetrOutsideShield(0));
4451     shieldindex=0;
4452   } else {
4453     Error("CorrectForShieldMaterial"," Wrong shield name\n");
4454     return 0;
4455   }
4456
4457   // do nothing if we are already beyond the shield
4458   Double_t rTrack = TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY());
4459   if(dir<0 && rTrack > rToGo) return 1; // going outward
4460   if(dir>0 && rTrack < rToGo) return 1; // going inward
4461
4462
4463   Double_t xToGo;
4464   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4465
4466   Int_t index=2*fCurrentEsdTrack+shieldindex;
4467
4468   Double_t xOverX0,x0,lengthTimesMeanDensity;
4469   Int_t nsteps=1;
4470
4471   switch(mode) {
4472   case 0:
4473     xOverX0 = AliITSRecoParam::Getdshield(shieldindex);
4474     x0 = AliITSRecoParam::GetX0shield(shieldindex);
4475     lengthTimesMeanDensity = xOverX0*x0;
4476     lengthTimesMeanDensity *= dir;
4477     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4478     break;
4479   case 1:
4480     nsteps= (Int_t)(TMath::Abs(t->GetX()-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4481     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4482     break;
4483   case 2:
4484     if(fxOverX0Shield[shieldindex]<0) BuildMaterialLUT("Shields");
4485     xOverX0 = fxOverX0Shield[shieldindex];
4486     lengthTimesMeanDensity = fxTimesRhoShield[shieldindex];
4487     lengthTimesMeanDensity *= dir;
4488     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4489     break;
4490   case 3:
4491     if(!fxOverX0ShieldTrks || index<0 || index>=2*fNtracks) Error("CorrectForShieldMaterial","Incorrect usage of UseTGeo option!\n");
4492     if(fxOverX0ShieldTrks[index]<0) {
4493       if (!t->PropagateToTGeo(xToGo,1,xOverX0,lengthTimesMeanDensity)) return 0;
4494       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4495                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4496       fxOverX0ShieldTrks[index] = TMath::Abs(xOverX0)/angle;
4497       fxTimesRhoShieldTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4498       return 1;
4499     }
4500     xOverX0 = fxOverX0ShieldTrks[index];
4501     lengthTimesMeanDensity = fxTimesRhoShieldTrks[index];
4502     lengthTimesMeanDensity *= dir;
4503     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4504     break;
4505   }
4506
4507   return 1;
4508 }
4509 //------------------------------------------------------------------------
4510 Int_t AliITStrackerMI::CorrectForLayerMaterial(AliITStrackMI *t,
4511                                                Int_t layerindex,
4512                                                Double_t oldGlobXYZ[3],
4513                                                TString direction) {
4514   //-------------------------------------------------------------------
4515   // Propagate beyond layer and correct for material
4516   // (material budget in different ways according to fUseTGeo value)
4517   // Add time if going outward (PropagateTo or PropagateToTGeo)
4518   //-------------------------------------------------------------------
4519
4520   // Define budget mode:
4521   // 0: material from AliITSRecoParam (hard coded)
4522   // 1: material from TGeo in stepsof X cm (on the fly)
4523   //    X = AliITSRecoParam::GetStepSizeTGeo()
4524   // 2: material from lut
4525   // 3: material from TGeo in one step (same for all hypotheses)
4526   Int_t mode;
4527   switch(fUseTGeo) {
4528   case 0:
4529     mode=0;
4530     break;
4531   case 1:
4532     mode=1;
4533     break;
4534   case 2:
4535     mode=2;
4536     break;
4537   case 3:
4538     if(fTrackingPhase.Contains("Clusters2Tracks"))
4539       { mode=3; } else { mode=1; }
4540     break;
4541   case 4:
4542     if(fTrackingPhase.Contains("Clusters2Tracks"))
4543       { mode=3; } else { mode=2; }
4544     break;
4545   default:
4546     mode=0;
4547     break;
4548   }
4549   if(fTrackingPhase.Contains("Default")) mode=0;
4550
4551   Float_t  dir = (direction.Contains("inward") ? 1. : -1.);
4552
4553   Double_t r=fgLayers[layerindex].GetR();
4554   Double_t deltar=(layerindex<2 ? 0.10*r : 0.05*r);
4555
4556   Double_t rToGo=TMath::Sqrt(t->GetX()*t->GetX()+t->GetY()*t->GetY())-deltar*dir;
4557   Double_t xToGo;
4558   if (!t->GetLocalXat(rToGo,xToGo)) return 0;
4559
4560   Int_t index=6*fCurrentEsdTrack+layerindex;
4561
4562
4563   Double_t xOverX0=0.0,x0=0.0,lengthTimesMeanDensity=0.0;
4564   Int_t nsteps=1;
4565
4566   // back before material (no correction)
4567   Double_t rOld,xOld;
4568   rOld=TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
4569   if (!t->GetLocalXat(rOld,xOld)) return 0;
4570   if (!t->Propagate(xOld)) return 0;
4571
4572   switch(mode) {
4573   case 0:
4574     xOverX0 = fgLayers[layerindex].GetThickness(t->GetY(),t->GetZ(),x0);
4575     lengthTimesMeanDensity = xOverX0*x0;
4576     lengthTimesMeanDensity *= dir;
4577     // Bring the track beyond the material
4578     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4579     break;
4580   case 1:
4581     nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4582     if (!t->PropagateToTGeo(xToGo,nsteps)) return 0; // cross the material and apply correction
4583     break;
4584   case 2:
4585     if(fxOverX0Layer[layerindex]<0) BuildMaterialLUT("Layers");
4586     xOverX0 = fxOverX0Layer[layerindex];
4587     lengthTimesMeanDensity = fxTimesRhoLayer[layerindex];
4588     lengthTimesMeanDensity *= dir;
4589     // Bring the track beyond the material
4590     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4591     break;
4592   case 3:
4593     if(!fxOverX0LayerTrks || index<0 || index>=6*fNtracks) Error("CorrectForLayerMaterial","Incorrect usage of UseTGeo option!\n");
4594     if(fxOverX0LayerTrks[index]<0) {
4595       nsteps = (Int_t)(TMath::Abs(xOld-xToGo)/AliITSReconstructor::GetRecoParam()->GetStepSizeTGeo())+1;
4596       if (!t->PropagateToTGeo(xToGo,nsteps,xOverX0,lengthTimesMeanDensity)) return 0;
4597       Double_t angle=TMath::Sqrt((1.+t->GetTgl()*t->GetTgl())/
4598                                  ((1.-t->GetSnp())*(1.+t->GetSnp())));
4599       fxOverX0LayerTrks[index] = TMath::Abs(xOverX0)/angle;
4600       fxTimesRhoLayerTrks[index] = TMath::Abs(lengthTimesMeanDensity)/angle;
4601       return 1;
4602     }
4603     xOverX0 = fxOverX0LayerTrks[index];
4604     lengthTimesMeanDensity = fxTimesRhoLayerTrks[index];
4605     lengthTimesMeanDensity *= dir;
4606     if (!t->PropagateTo(xToGo,xOverX0,lengthTimesMeanDensity/xOverX0)) return 0;
4607     break;
4608   }
4609
4610
4611   return 1;
4612 }
4613 //------------------------------------------------------------------------
4614 void AliITStrackerMI::MakeTrksMaterialLUT(Int_t ntracks) {
4615   //-----------------------------------------------------------------
4616   // Initialize LUT for storing material for each prolonged track
4617   //-----------------------------------------------------------------
4618   fxOverX0PipeTrks = new Float_t[ntracks];
4619   fxTimesRhoPipeTrks = new Float_t[ntracks];
4620   fxOverX0ShieldTrks = new Float_t[ntracks*2];
4621   fxTimesRhoShieldTrks = new Float_t[ntracks*2];
4622   fxOverX0LayerTrks = new Float_t[ntracks*6];
4623   fxTimesRhoLayerTrks = new Float_t[ntracks*6];
4624
4625   for(Int_t i=0; i<ntracks; i++) {
4626     fxOverX0PipeTrks[i] = -1.;
4627     fxTimesRhoPipeTrks[i] = -1.;
4628   }
4629   for(Int_t j=0; j<ntracks*2; j++) {
4630     fxOverX0ShieldTrks[j] = -1.;
4631     fxTimesRhoShieldTrks[j] = -1.;
4632   }
4633   for(Int_t k=0; k<ntracks*6; k++) {
4634     fxOverX0LayerTrks[k] = -1.;
4635     fxTimesRhoLayerTrks[k] = -1.;
4636   }
4637
4638   fNtracks = ntracks;
4639
4640   return;
4641 }
4642 //------------------------------------------------------------------------
4643 void AliITStrackerMI::DeleteTrksMaterialLUT() {
4644   //-----------------------------------------------------------------
4645   // Delete LUT for storing material for each prolonged track
4646   //-----------------------------------------------------------------
4647   if(fxOverX0PipeTrks) {
4648     delete [] fxOverX0PipeTrks; fxOverX0PipeTrks = 0;
4649   }
4650   if(fxOverX0ShieldTrks) {
4651     delete [] fxOverX0ShieldTrks; fxOverX0ShieldTrks = 0;
4652   }
4653
4654   if(fxOverX0LayerTrks) {
4655     delete [] fxOverX0LayerTrks;  fxOverX0LayerTrks = 0;
4656   }
4657   if(fxTimesRhoPipeTrks) {
4658     delete [] fxTimesRhoPipeTrks;  fxTimesRhoPipeTrks = 0;
4659   }
4660   if(fxTimesRhoShieldTrks) {
4661     delete [] fxTimesRhoShieldTrks; fxTimesRhoShieldTrks = 0;
4662   }
4663   if(fxTimesRhoLayerTrks) {
4664     delete [] fxTimesRhoLayerTrks; fxTimesRhoLayerTrks = 0;
4665   }
4666   return;
4667 }
4668 //------------------------------------------------------------------------
4669 void AliITStrackerMI::SetForceSkippingOfLayer() {
4670   //-----------------------------------------------------------------
4671   // Check if we are forced to skip layers
4672   // either we set to skip them in RecoParam
4673   // or they were off during data-taking
4674   //-----------------------------------------------------------------
4675
4676   const AliEventInfo *eventInfo = GetEventInfo();
4677
4678   for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4679     fForceSkippingOfLayer[l] = 0;
4680     // check reco param
4681     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4682     // check run info
4683
4684     if(eventInfo &&
4685        AliITSReconstructor::GetRecoParam()->GetSkipSubdetsNotInTriggerCluster()) {
4686       AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4687       if(l==0 || l==1)  {
4688         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4689       } else if(l==2 || l==3) {
4690         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4691       } else {
4692         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4693       }
4694     }
4695   }
4696   return;
4697 }
4698 //------------------------------------------------------------------------
4699 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4700                                       Int_t ilayer,Int_t idet) const {
4701   //-----------------------------------------------------------------
4702   // This method is used to decide whether to allow a prolongation
4703   // without clusters, because we want to skip the layer.
4704   // In this case the return value is > 0:
4705   // return 1: the user requested to skip a layer
4706   // return 2: track outside z acceptance
4707   //-----------------------------------------------------------------
4708
4709   if (ForceSkippingOfLayer(ilayer)) return 1;
4710
4711   Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4712
4713   if (idet<0 &&  // out in z
4714       ilayer>innerLayCanSkip &&
4715       AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4716     // check if track will cross SPD outer layer
4717     Double_t phiAtSPD2,zAtSPD2;
4718     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4719       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4720     }
4721     return 2; // always allow skipping, changed on 05.11.2009
4722   }
4723
4724   return 0;
4725 }
4726 //------------------------------------------------------------------------
4727 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4728                                      Int_t ilayer,Int_t idet,
4729                                      Double_t dz,Double_t dy,
4730                                      Bool_t noClusters) const {
4731   //-----------------------------------------------------------------
4732   // This method is used to decide whether to allow a prolongation
4733   // without clusters, because there is a dead zone in the road.
4734   // In this case the return value is > 0:
4735   // return 1: dead zone at z=0,+-7cm in SPD
4736   //     This method assumes that fSPDdetzcentre is ordered from -z to +z
4737   // return 2: all road is "bad" (dead or noisy) from the OCDB
4738   // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4739   // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4740   //-----------------------------------------------------------------
4741
4742   // check dead zones at z=0,+-7cm in the SPD
4743   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4744     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4745                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4746                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4747     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4748                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4749                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4750     for (Int_t i=0; i<3; i++)
4751       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4752         AliDebug(2,Form("crack SPD %d track z %f   %f   %f  %f\n",ilayer,track->GetZ(),dz,zmaxdead[i],zmindead[i]));
4753         if (GetSPDDeadZoneProbability(track->GetZ(),TMath::Sqrt(track->GetSigmaZ2()))>0.1) return 1;
4754       }
4755   }
4756
4757   // check bad zones from OCDB
4758   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4759
4760   if (idet<0) return 0;
4761
4762   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4763
4764   Int_t detType=-1;
4765   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4766   if (ilayer==0 || ilayer==1) {        // ----------  SPD
4767     detType = 0;
4768   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
4769     detType = 1;
4770     detSizeFactorX *= 2.;
4771   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
4772     detType = 2;
4773   }
4774   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4775   if (detType==2) segm->SetLayer(ilayer+1);
4776   Float_t detSizeX = detSizeFactorX*segm->Dx();
4777   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4778
4779   // check if the road overlaps with bad chips
4780   Float_t xloc,zloc;
4781   if(!(LocalModuleCoord(ilayer,idet,track,xloc,zloc)))return 0;
4782   Float_t zlocmin = zloc-dz;
4783   Float_t zlocmax = zloc+dz;
4784   Float_t xlocmin = xloc-dy;
4785   Float_t xlocmax = xloc+dy;
4786   Int_t chipsInRoad[100];
4787
4788   // check if road goes out of detector
4789   Bool_t touchNeighbourDet=kFALSE;
4790   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4791   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.4999*detSizeX; touchNeighbourDet=kTRUE;}
4792   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4793   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.4999*detSizeZ; touchNeighbourDet=kTRUE;}
4794   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));
4795
4796   // check if this detector is bad
4797   if (det.IsBad()) {
4798     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
4799     if(!touchNeighbourDet) {
4800       return 2; // all detectors in road are bad
4801     } else {
4802       return 3; // at least one is bad
4803     }
4804   }
4805
4806   if(zlocmin>zlocmax)return 0;
4807   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4808   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4809   if (!nChipsInRoad) return 0;
4810
4811   Bool_t anyBad=kFALSE,anyGood=kFALSE;
4812   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4813     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4814     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4815     if (det.IsChipBad(chipsInRoad[iCh])) {
4816       anyBad=kTRUE;
4817     } else {
4818       anyGood=kTRUE;
4819     }
4820   }
4821
4822   if (!anyGood) {
4823     if(!touchNeighbourDet) {
4824       AliDebug(2,"all bad in road");
4825       return 2;  // all chips in road are bad
4826     } else {
4827       return 3; // at least a bad chip in road
4828     }
4829   }
4830
4831   if (anyBad) {
4832     AliDebug(2,"at least a bad in road");
4833     return 3; // at least a bad chip in road
4834   }
4835
4836
4837   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4838       || !noClusters) return 0;
4839
4840   // There are no clusters in road: check if there is at least
4841   // a bad SPD pixel or SDD anode or SSD strips on both sides
4842
4843   Int_t idetInITS=idet;
4844   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4845
4846   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4847     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4848     return 4;
4849   }
4850   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4851
4852   return 0;
4853 }
4854 //------------------------------------------------------------------------
4855 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4856                                        const AliITStrackMI *track,
4857                                        Float_t &xloc,Float_t &zloc) const {
4858   //-----------------------------------------------------------------
4859   // Gives position of track in local module ref. frame
4860   //-----------------------------------------------------------------
4861
4862   xloc=0.;
4863   zloc=0.;
4864
4865   if(idet<0) return kTRUE; // track out of z acceptance of layer
4866
4867   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4868
4869   Int_t lad = Int_t(idet/ndet) + 1;
4870
4871   Int_t det = idet - (lad-1)*ndet + 1;
4872
4873   Double_t xyzGlob[3],xyzLoc[3];
4874
4875   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4876   // take into account the misalignment: xyz at real detector plane
4877   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4878
4879   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4880
4881   xloc = (Float_t)xyzLoc[0];
4882   zloc = (Float_t)xyzLoc[2];
4883
4884   return kTRUE;
4885 }
4886 //------------------------------------------------------------------------
4887 //------------------------------------------------------------------------
4888 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer){
4889 //
4890 // Method to be optimized further:
4891 // Aim: decide whether a track can be used for PlaneEff evaluation
4892 //      the decision is taken based on the track quality at the layer under study
4893 //      no information on the clusters on this layer has to be used
4894 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4895 //      the cut is done on number of sigmas from the boundaries
4896 //
4897 //  Input: Actual track, layer [0,5] under study
4898 //  Output: none
4899 //  Return: kTRUE if this is a good track
4900 //
4901 // it will apply a pre-selection to obtain good quality tracks.
4902 // Here also  you will have the possibility to put a control on the
4903 // impact point of the track on the basic block, in order to exclude border regions
4904 // this will be done by calling a proper method of the AliITSPlaneEff class.
4905 //
4906 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4907 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
4908 //
4909   Int_t index[AliITSgeomTGeo::kNLayers];
4910   Int_t k;
4911   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4912   //
4913   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4914     index[k]=clusters[k];
4915   }
4916
4917   if(!fPlaneEff)
4918     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4919   AliITSlayer &layer=fgLayers[ilayer];
4920   Double_t r=layer.GetR();
4921   AliITStrackMI tmp(*track);
4922
4923 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4924   Int_t nclout=0; Int_t nclin=0;
4925   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) { // count n. of cluster in outermost layers
4926  AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4927                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4928    // if (tmp.GetClIndex(lay)>=0) nclout++;
4929 if(index[lay]>=0)nclout++;
4930   }
4931   for(Int_t lay=ilayer-1; lay>=0;lay--) { // count n. of cluster in innermost layers
4932     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4933                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4934    if (index[lay]>=0) nclin++;
4935   }
4936   Int_t ncl=nclout+nclin;
4937   Bool_t nextout = kFALSE;
4938   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4939   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4940   Bool_t nextin = kFALSE;
4941   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4942   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4943   // maximum number of missing clusters allowed in outermost layers
4944   if(nclout<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersOutPlaneEff())
4945      return kFALSE;
4946   // maximum number of missing clusters allowed (both in innermost and in outermost layers)
4947   if(ncl<AliITSgeomTGeo::kNLayers-1-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4948      return kFALSE;
4949   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
4950   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
4951   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4952  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
4953
4954 // detector number
4955   Double_t phi,z;
4956   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4957   Int_t idet=layer.FindDetectorIndex(phi,z);
4958   if(idet<0) { AliInfo(Form("cannot find detector"));
4959     return kFALSE;}
4960
4961   // here check if it has good Chi Square.
4962
4963   //propagate to the intersection with the detector plane
4964   const AliITSdetector &det=layer.GetDetector(idet);
4965   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4966
4967   Float_t locx; //
4968   Float_t locz; //
4969   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4970   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4971   if(key>fPlaneEff->Nblock()) return kFALSE;
4972   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4973   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4974   //***************
4975   // DEFINITION OF SEARCH ROAD FOR accepting a track
4976   //
4977   Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXFromBoundaryPlaneEff();
4978   Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZFromBoundaryPlaneEff();
4979   Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXFromBoundaryPlaneEff();
4980   Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZFromBoundaryPlaneEff();
4981   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
4982   // those are precisions in the tracking reference system
4983   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
4984   // Use it also for the module reference system, as it is done for RecPoints
4985
4986   if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigFrmBndPlaneEff()){
4987    if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) dx -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
4988    else dx -= distx;
4989
4990    if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) dz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
4991    else dz -= distz;
4992   }
4993
4994   // exclude tracks at boundary between detectors
4995   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4996   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4997   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4998   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
4999   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
5000   if ( (locx-dx < blockXmn+boundaryWidth) ||
5001        (locx+dx > blockXmx-boundaryWidth) ||
5002        (locz-dz < blockZmn+boundaryWidth) ||
5003        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
5004
5005   if(ilayer==0){
5006        const AliESDEvent *myesd = ((AliESDtrack*)tmp.GetESDtrack())->GetESDEvent();
5007        //The beam pipe
5008        if (CorrectForPipeMaterial(&tmp,"inward")) {
5009           const AliESDVertex* vtx = myesd->GetVertex();
5010           if(!vtx) return kFALSE;
5011           Double_t ddz[2],cov[3];
5012           Double_t maxD=3.;
5013           if(!tmp.PropagateToDCA(vtx,AliTracker::GetBz(),maxD,ddz,cov)) return kFALSE;
5014           if(TMath::Abs(ddz[0])>=AliITSReconstructor::GetRecoParam()->GetDCACutPlaneEff()) return kFALSE;
5015
5016           Double_t covar[6]; vtx->GetCovMatrix(covar);
5017           Double_t p[2]={tmp.GetParameter()[0]-ddz[0],tmp.GetParameter()[1]-ddz[1]};
5018           Double_t c[3]={covar[2],0.,covar[5]};
5019           Double_t chi2= ((AliESDtrack*)tmp.GetESDtrack())->GetPredictedChi2(p,c);
5020           if (chi2>AliITSReconstructor::GetRecoParam()->GetVertexChi2CutPlaneEff()) return kFALSE;  // Use this to cut on chi^2
5021
5022        } else return kFALSE;
5023     }
5024
5025
5026   return kTRUE;
5027 }
5028 //------------------------------------------------------------------------
5029 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
5030 //
5031 // This Method has to be optimized! For the time-being it uses the same criteria
5032 // as those used in the search of extra clusters for overlapping modules.
5033 //
5034 // Method Purpose: estabilish whether a track has produced a recpoint or not
5035 //                 in the layer under study (For Plane efficiency)
5036 //
5037 // inputs: AliITStrackMI* track  (pointer to a usable track)
5038 // outputs: none
5039 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
5040 //               traversed by this very track. In details:
5041 //               - if a cluster can be associated to the track then call
5042 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
5043 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
5044 //
5045   if(!fPlaneEff)
5046     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
5047   AliITSlayer &layer=fgLayers[ilayer];
5048   Double_t r=layer.GetR();
5049   AliITStrackMI tmp(*track);
5050
5051 // detector number
5052   Double_t phi,z;
5053   if (!tmp.GetPhiZat(r,phi,z)) return;
5054   Int_t idet=layer.FindDetectorIndex(phi,z);
5055
5056   if(idet<0) { AliInfo(Form("cannot find detector"));
5057     return;}
5058
5059
5060 //propagate to the intersection with the detector plane
5061   const AliITSdetector &det=layer.GetDetector(idet);
5062   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
5063
5064
5065 //***************
5066 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
5067 //
5068   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
5069                     TMath::Sqrt(tmp.GetSigmaZ2() +
5070                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5071                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5072                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
5073   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
5074                     TMath::Sqrt(tmp.GetSigmaY2() +
5075                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5076                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5077                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
5078
5079 // road in global (rphi,z) [i.e. in tracking ref. system]
5080   Double_t zmin = tmp.GetZ() - dz;
5081   Double_t zmax = tmp.GetZ() + dz;
5082   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
5083   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
5084
5085 // select clusters in road
5086   layer.SelectClusters(zmin,zmax,ymin,ymax);
5087
5088 // Define criteria for track-cluster association
5089   Double_t msz = tmp.GetSigmaZ2() +
5090   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5091   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
5092   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
5093   Double_t msy = tmp.GetSigmaY2() +
5094   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5095   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
5096   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
5097   if (tmp.GetConstrain()) {
5098     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
5099     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
5100   }  else {
5101     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
5102     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
5103   }
5104
5105   if(AliITSReconstructor::GetRecoParam()->GetSwitchOffStdSearchClusPlaneEff()){
5106      Double_t nsigx=AliITSReconstructor::GetRecoParam()->GetNSigXSearchClusterPlaneEff();
5107      Double_t nsigz=AliITSReconstructor::GetRecoParam()->GetNSigZSearchClusterPlaneEff();
5108      Double_t distx=AliITSReconstructor::GetRecoParam()->GetDistXSearchClusterPlaneEff();
5109      Double_t distz=AliITSReconstructor::GetRecoParam()->GetDistZSearchClusterPlaneEff();
5110      msy = nsigx*TMath::Sqrt(tmp.GetSigmaY2()) + distx;
5111      msz = nsigz*TMath::Sqrt(tmp.GetSigmaZ2()) + distz;
5112
5113   if(AliITSReconstructor::GetRecoParam()->GetSwitchOnMaxDistNSigSrhClusPlaneEff()){
5114      if(nsigx*TMath::Sqrt(tmp.GetSigmaY2())<=distx) msy -= nsigx*TMath::Sqrt(tmp.GetSigmaY2());
5115      else msy -= distx;
5116
5117      if(nsigz*TMath::Sqrt(tmp.GetSigmaZ2())<=distz) msz -= nsigz*TMath::Sqrt(tmp.GetSigmaZ2());
5118      else msz -= distz;
5119   }
5120
5121   msy *= msy;
5122   msz *= msz;
5123
5124   }
5125
5126   if(msz==0 || msy==0){AliWarning("UseTrackForPlaneEff: null search frame"); return;}
5127
5128   msz = 1./msz; // 1/RoadZ^2
5129   msy = 1./msy; // 1/RoadY^2
5130
5131   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
5132   Int_t idetc=-1;
5133   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
5134   //Double_t  tolerance=0.2;
5135   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
5136     idetc = cl->GetDetectorIndex();
5137     if(idet!=idetc) continue;
5138     //Int_t ilay = cl->GetLayer();
5139
5140     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
5141     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
5142
5143     Double_t chi2=tmp.GetPredictedChi2(cl);
5144     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
5145   }*/
5146   Float_t locx; //
5147   Float_t locz; //
5148   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
5149 //
5150   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
5151   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
5152   if(key>fPlaneEff->Nblock()) return;
5153   Bool_t found=kFALSE;
5154   //if (ci>=0) {
5155   Double_t chi2;
5156   while ((cl=layer.GetNextCluster(clidx))!=0) {
5157     idetc = cl->GetDetectorIndex();
5158     if(idet!=idetc) continue;
5159     // here real control to see whether the cluster can be associated to the track.
5160     // cluster not associated to track
5161     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
5162          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
5163     // calculate track-clusters chi2
5164     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
5165                                                // in particular, the error associated to the cluster
5166     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
5167     // chi2 cut
5168     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
5169     found=kTRUE;
5170     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
5171    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
5172    // track->SetExtraModule(ilayer,idetExtra);
5173   }
5174   if(!fPlaneEff->UpDatePlaneEff(found,key))
5175        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
5176
5177 // this for FO efficiency studies (only for SPD) //
5178    UInt_t keyFO=999999;
5179    Bool_t foundFO=kFALSE;
5180    if(ilayer<2){ //ONLY SPD layers for FastOr studies
5181     TBits mapFO = fkDetTypeRec->GetFastOrFiredMap();
5182     Int_t phase = (fEsd->GetBunchCrossNumber())%4;
5183     if(!fSPDChipIntPlaneEff[key]){
5184       AliITSPlaneEffSPD spd;
5185       keyFO = spd.SwitchChipKeyNumbering(key);
5186       if(mapFO.TestBitNumber(keyFO))foundFO=kTRUE;
5187        keyFO = key + (AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip)*(phase+1);
5188        if(keyFO<AliITSPlaneEffSPD::kNModule*AliITSPlaneEffSPD::kNChip) {
5189          AliWarning(Form("UseTrackForPlaneEff: too small keyF0 (= %d), setting it to 999999",keyFO));
5190          keyFO=999999;
5191        }
5192        if(!fPlaneEff->UpDatePlaneEff(foundFO,keyFO))
5193           AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for FastOR for key=%d",keyFO));
5194      }
5195   }
5196
5197
5198
5199   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
5200     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
5201     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
5202     Int_t cltype[2]={-999,-999};
5203                                                           // and the module
5204
5205     Float_t angleModTrack[3]={99999.,99999.,99999.}; // angles (phi, z and "absolute angle") between the track and the mormal to the module (see below)
5206
5207     tr[0]=locx;
5208     tr[1]=locz;
5209     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
5210     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
5211
5212     if (found){
5213       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
5214       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
5215       cltype[0]=layer.GetCluster(ci)->GetNy();
5216       cltype[1]=layer.GetCluster(ci)->GetNz();
5217
5218      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
5219      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
5220      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
5221      // It is computed properly by calling the method
5222      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
5223      // T
5224      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
5225       //if (tmp.PropagateTo(x,0.,0.)) {
5226         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
5227         AliCluster c(*layer.GetCluster(ci));
5228         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
5229         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
5230         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
5231         clu[2]=TMath::Sqrt(c.GetSigmaY2());
5232         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
5233       //}
5234     }
5235   // Compute the angles between the track and the module
5236       // compute the angle "in phi direction", i.e. the angle in the transverse plane
5237       // between the normal to the module and the projection (in the transverse plane) of the
5238       // track trajectory
5239     // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
5240     Float_t tgl = tmp.GetTgl();
5241     Float_t phitr   = tmp.GetSnp();
5242     phitr = TMath::ASin(phitr);
5243     Int_t volId = AliGeomManager::LayerToVolUIDSafe(ilayer+1 ,idet );
5244
5245     Double_t tra[3]; AliGeomManager::GetOrigTranslation(volId,tra);
5246     Double_t rot[9]; AliGeomManager::GetOrigRotation(volId,rot);
5247    Double_t alpha =0.;
5248     alpha = tmp.GetAlpha();
5249     Double_t phiglob = alpha+phitr;
5250     Double_t p[3];
5251     p[0] = TMath::Cos(phiglob);
5252     p[1] = TMath::Sin(phiglob);
5253     p[2] = tgl;
5254     TVector3 pvec(p[0],p[1],p[2]);
5255     TVector3 normvec(rot[1],rot[4],rot[7]);
5256     Double_t angle = pvec.Angle(normvec);
5257
5258     if(angle>0.5*TMath::Pi()) angle = (TMath::Pi()-angle);
5259     angle *= 180./TMath::Pi();
5260
5261     //Trasverse Plane
5262     TVector3 pt(p[0],p[1],0);
5263     TVector3 normt(rot[1],rot[4],0);
5264     Double_t anglet = pt.Angle(normt);
5265
5266     Double_t phiPt = TMath::ATan2(p[1],p[0]);
5267     if(phiPt<0)phiPt+=2.*TMath::Pi();
5268     Double_t phiNorm = TMath::ATan2(rot[4],rot[1]);
5269     if(phiNorm<0) phiNorm+=2.*TMath::Pi();
5270     if(anglet>0.5*TMath::Pi()) anglet = (TMath::Pi()-anglet);
5271     if(phiNorm>phiPt) anglet*=-1.;// pt-->normt  clockwise: anglet>0
5272     if((phiNorm-phiPt)>TMath::Pi()) anglet*=-1.;
5273     anglet *= 180./TMath::Pi();
5274
5275      angleModTrack[2]=(Float_t) angle;
5276      angleModTrack[0]=(Float_t) anglet;
5277      // now the "angle in z" (much easier, i.e. the angle between the z axis and the track momentum + 90)
5278     angleModTrack[1]=TMath::ACos(tgl/TMath::Sqrt(tgl*tgl+1.));
5279     angleModTrack[1]-=TMath::Pi()/2.; // range of angle is -pi/2 , pi/2
5280     angleModTrack[1]*=180./TMath::Pi(); // in degree
5281
5282     fPlaneEff->FillHistos(key,found,tr,clu,cltype,angleModTrack);
5283
5284     // For FO efficiency studies of SPD
5285     if(ilayer<2 && !fSPDChipIntPlaneEff[key]) fPlaneEff->FillHistos(keyFO,foundFO,tr,clu,cltype,angleModTrack);
5286   }
5287   if(ilayer<2) fSPDChipIntPlaneEff[key]=kTRUE;
5288 return;
5289 }
5290
5291 Int_t AliITStrackerMI::FindClusterOfTrack(int label, int lr, int* store) const //RS
5292 {
5293   // find the MC cluster for the label
5294   return fgLayers[lr].FindClusterForLabel(label,store);
5295 }
5296
5297 /*
5298 Int_t AliITStrackerMI::GetPattern(const AliITStrackMI* track, char* patt)
5299 {
5300   // creates pattarn of hits marking fake/corrects by f/c. Used for debugging (RS)
5301   strncpy(patt,"......",6);
5302   int tpcLabel = 0;
5303   if (track->GetESDtrack()) tpcLabel = track->GetESDtrack()->GetTPCLabel();
5304   //
5305   int nwrong = 0;
5306   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
5307     Int_t cindex = track->GetClusterIndex(i);
5308     Int_t l=(cindex & 0xf0000000) >> 28;
5309     AliITSRecPoint *cl = (AliITSRecPoint*)GetCluster(cindex);
5310     Int_t isWrong=1;
5311     for (Int_t ind=0;ind<3;ind++) if (cl->GetLabel(ind)==TMath::Abs(tpcLabel)) isWrong=0;
5312     patt[l] = isWrong ? 'f':'c';
5313     nwrong+=isWrong;
5314   }
5315   return nwrong;
5316 }
5317 */
5318 //------------------------------------------------------------------------
5319 Int_t AliITStrackerMI::AliITSlayer::FindClusterForLabel(Int_t label, Int_t *store) const
5320 { //RS
5321   //--------------------------------------------------------------------
5322
5323   int nfound = 0;
5324   for (int ic=0;ic<fN;ic++) {
5325     const AliITSRecPoint *cl = GetCluster(ic);
5326     for (int i=0;i<3;i++) if (cl->GetLabel(i)==label) {
5327         if (nfound<50) {
5328           if (store) store[nfound] = ic;
5329           nfound++;
5330         }
5331         break;
5332       }
5333   }
5334   return nfound;
5335 }