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