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