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