ITS/AliITStrackerMI.cxx

   1 /**************************************************************************
   2  * Copyright(c) 2007-2009, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id$ */
  16
  17 //-------------------------------------------------------------------------
  18 //               Implementation of the ITS tracker class
  19 //    It reads AliITSRecPoint clusters and creates AliITStrackMI tracks
  20 //                   and fills with them the ESD
  21 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  22 //          Current support and development:
  23 //                     Andrea Dainese, andrea.dainese@lnl.infn.it
  24 //     dE/dx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  25 //     Params moved to AliITSRecoParam by: Andrea Dainese, INFN
  26 //     Material budget from TGeo by: Ludovic Gaudichet & Andrea Dainese, INFN
  27 //-------------------------------------------------------------------------
  28
  29 #include <TMatrixD.h>
  30 #include <TTree.h>
  31 #include <TDatabasePDG.h>
  32 #include <TString.h>
  33 #include <TRandom.h>
  34 #include <TTreeStream.h>
  35
  36 #include "AliLog.h"
  37 #include "AliITSPlaneEff.h"
  38 #include "AliITSCalibrationSPD.h"
  39 #include "AliITSCalibrationSDD.h"
  40 #include "AliITSCalibrationSSD.h"
  41 #include "AliCDBEntry.h"
  42 #include "AliCDBManager.h"
  43 #include "AliAlignObj.h"
  44 #include "AliTrackPointArray.h"
  45 #include "AliESDVertex.h"
  46 #include "AliESDEvent.h"
  47 #include "AliESDtrack.h"
  48 #include "AliV0.h"
  49 #include "AliHelix.h"
  50 #include "AliITSChannelStatus.h"
  51 #include "AliITSDetTypeRec.h"
  52 #include "AliITSRecPoint.h"
  53 #include "AliITSgeomTGeo.h"
  54 #include "AliITSReconstructor.h"
  55 #include "AliITSClusterParam.h"
  56 #include "AliITSsegmentation.h"
  57 #include "AliITSCalibration.h"
  58 #include "AliITSPlaneEffSPD.h"
  59 #include "AliITSPlaneEffSDD.h"
  60 #include "AliITSPlaneEffSSD.h"
  61 #include "AliITSV0Finder.h"
  62 #include "AliITStrackerMI.h"
  63
  64 ClassImp(AliITStrackerMI)
  65
  66 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[AliITSgeomTGeo::kNLayers]; // ITS layers
  67
  68 AliITStrackerMI::AliITStrackerMI():AliTracker(),
  69 fI(0),
  70 fBestTrack(),
  71 fTrackToFollow(),
  72 fTrackHypothesys(),
  73 fBestHypothesys(),
  74 fOriginal(),
  75 fCurrentEsdTrack(),
  76 fPass(0),
  77 fAfterV0(kFALSE),
  78 fLastLayerToTrackTo(0),
  79 fCoefficients(0),
  80 fEsd(0),
  81 fTrackingPhase("Default"),
  82 fUseTGeo(3),
  83 fNtracks(0),
  84 fxOverX0Pipe(-1.),
  85 fxTimesRhoPipe(-1.),
  86 fxOverX0PipeTrks(0),
  87 fxTimesRhoPipeTrks(0),
  88 fxOverX0ShieldTrks(0),
  89 fxTimesRhoShieldTrks(0),
  90 fxOverX0LayerTrks(0),
  91 fxTimesRhoLayerTrks(0),
  92 fDebugStreamer(0),
  93 fITSChannelStatus(0),
  94 fkDetTypeRec(0),
  95 fPlaneEff(0) {
  96   //Default constructor
  97   Int_t i;
  98   for(i=0;i<4;i++) fSPDdetzcentre[i]=0.;
  99   for(i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 100   for(i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 101 }
 102 //------------------------------------------------------------------------
 103 AliITStrackerMI::AliITStrackerMI(const Char_t *geom) : AliTracker(),
 104 fI(AliITSgeomTGeo::GetNLayers()),
 105 fBestTrack(),
 106 fTrackToFollow(),
 107 fTrackHypothesys(),
 108 fBestHypothesys(),
 109 fOriginal(),
 110 fCurrentEsdTrack(),
 111 fPass(0),
 112 fAfterV0(kFALSE),
 113 fLastLayerToTrackTo(AliITSRecoParam::GetLastLayerToTrackTo()),
 114 fCoefficients(0),
 115 fEsd(0),
 116 fTrackingPhase("Default"),
 117 fUseTGeo(3),
 118 fNtracks(0),
 119 fxOverX0Pipe(-1.),
 120 fxTimesRhoPipe(-1.),
 121 fxOverX0PipeTrks(0),
 122 fxTimesRhoPipeTrks(0),
 123 fxOverX0ShieldTrks(0),
 124 fxTimesRhoShieldTrks(0),
 125 fxOverX0LayerTrks(0),
 126 fxTimesRhoLayerTrks(0),
 127 fDebugStreamer(0),
 128 fITSChannelStatus(0),
 129 fkDetTypeRec(0),
 130 fPlaneEff(0) {
 131   //--------------------------------------------------------------------
 132   //This is the AliITStrackerMI constructor
 133   //--------------------------------------------------------------------
 134   if (geom) {
 135     AliWarning("\"geom\" is actually a dummy argument !");
 136   }
 137
 138   fCoefficients = 0;
 139   fAfterV0     = kFALSE;
 140
 141   for (Int_t i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 142     Int_t nlad=AliITSgeomTGeo::GetNLadders(i);
 143     Int_t ndet=AliITSgeomTGeo::GetNDetectors(i);
 144
 145     Double_t xyz[3], &x=xyz[0], &y=xyz[1], &z=xyz[2];
 146     AliITSgeomTGeo::GetOrigTranslation(i,1,1,xyz);
 147     Double_t poff=TMath::ATan2(y,x);
 148     Double_t zoff=z;
 149     Double_t r=TMath::Sqrt(x*x + y*y);
 150
 151     AliITSgeomTGeo::GetOrigTranslation(i,1,2,xyz);
 152     r += TMath::Sqrt(x*x + y*y);
 153     AliITSgeomTGeo::GetOrigTranslation(i,2,1,xyz);
 154     r += TMath::Sqrt(x*x + y*y);
 155     AliITSgeomTGeo::GetOrigTranslation(i,2,2,xyz);
 156     r += TMath::Sqrt(x*x + y*y);
 157     r*=0.25;
 158
 159     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
 160
 161     for (Int_t j=1; j<nlad+1; j++) {
 162       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
 163         TGeoHMatrix m; AliITSgeomTGeo::GetOrigMatrix(i,j,k,m);
 164         const TGeoHMatrix *tm=AliITSgeomTGeo::GetTracking2LocalMatrix(i,j,k);
 165         m.Multiply(tm);
 166         Double_t txyz[3]={0.};
 167         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 168         m.LocalToMaster(txyz,xyz);
 169         r=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 170         Double_t phi=TMath::ATan2(xyz[1],xyz[0]);
 171
 172         if (phi<0) phi+=TMath::TwoPi();
 173         else if (phi>=TMath::TwoPi()) phi-=TMath::TwoPi();
 174
 175         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 176         new(&det) AliITSdetector(r,phi);
 177         // compute the real radius (with misalignment)
 178         TGeoHMatrix mmisal(*(AliITSgeomTGeo::GetMatrix(i,j,k)));
 179         mmisal.Multiply(tm);
 180         xyz[0]=0.;xyz[1]=0.;xyz[2]=0.;
 181         mmisal.LocalToMaster(txyz,xyz);
 182         Double_t rmisal=TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]*xyz[1]);
 183         det.SetRmisal(rmisal);
 184
 185       } // end loop on detectors
 186     } // end loop on ladders
 187     fForceSkippingOfLayer[i] = 0;
 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   fLastLayerToTrackTo=AliITSRecoParam::GetLastLayerToTrackTo();
 205   for (Int_t i=0;i<100000;i++){
 206     fBestTrackIndex[i]=0;
 207   }
 208
 209   // store positions of centre of SPD modules (in z)
 210   Double_t tr[3];
 211   AliITSgeomTGeo::GetTranslation(1,1,1,tr);
 212   fSPDdetzcentre[0] = tr[2];
 213   AliITSgeomTGeo::GetTranslation(1,1,2,tr);
 214   fSPDdetzcentre[1] = tr[2];
 215   AliITSgeomTGeo::GetTranslation(1,1,3,tr);
 216   fSPDdetzcentre[2] = tr[2];
 217   AliITSgeomTGeo::GetTranslation(1,1,4,tr);
 218   fSPDdetzcentre[3] = tr[2];
 219
 220   fUseTGeo = AliITSReconstructor::GetRecoParam()->GetUseTGeoInTracker();
 221   if(AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance() && fUseTGeo!=1 && fUseTGeo!=3) {
 222     AliWarning("fUseTGeo changed to 3 because fExtendedEtaAcceptance is kTRUE");
 223     fUseTGeo = 3;
 224   }
 225
 226   for(Int_t i=0;i<2;i++) {fxOverX0Shield[i]=-1.;fxTimesRhoShield[i]=-1.;}
 227   for(Int_t i=0;i<6;i++) {fxOverX0Layer[i]=-1.;fxTimesRhoLayer[i]=-1.;}
 228
 229   fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
 230
 231   // only for plane efficiency evaluation
 232   if (AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 233       AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0) {
 234     Int_t iplane=AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff();
 235     if(!AliITSReconstructor::GetRecoParam()->GetLayersToSkip(iplane)==1)
 236       AliWarning(Form("Evaluation of Plane Eff for layer %d will be attempted without removing it from tracker",iplane));
 237     if (iplane<2) fPlaneEff = new AliITSPlaneEffSPD();
 238     else if (iplane<4) fPlaneEff = new AliITSPlaneEffSDD();
 239     else fPlaneEff = new AliITSPlaneEffSSD();
 240     if(AliITSReconstructor::GetRecoParam()->GetReadPlaneEffFromOCDB())
 241        if(!fPlaneEff->ReadFromCDB()) {AliWarning("AliITStrackerMI reading of AliITSPlaneEff from OCDB failed") ;}
 242     if(AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) fPlaneEff->SetCreateHistos(kTRUE);
 243   }
 244 }
 245 //------------------------------------------------------------------------
 246 AliITStrackerMI::AliITStrackerMI(const AliITStrackerMI &tracker):AliTracker(tracker),
 247 fI(tracker.fI),
 248 fBestTrack(tracker.fBestTrack),
 249 fTrackToFollow(tracker.fTrackToFollow),
 250 fTrackHypothesys(tracker.fTrackHypothesys),
 251 fBestHypothesys(tracker.fBestHypothesys),
 252 fOriginal(tracker.fOriginal),
 253 fCurrentEsdTrack(tracker.fCurrentEsdTrack),
 254 fPass(tracker.fPass),
 255 fAfterV0(tracker.fAfterV0),
 256 fLastLayerToTrackTo(tracker.fLastLayerToTrackTo),
 257 fCoefficients(tracker.fCoefficients),
 258 fEsd(tracker.fEsd),
 259 fTrackingPhase(tracker.fTrackingPhase),
 260 fUseTGeo(tracker.fUseTGeo),
 261 fNtracks(tracker.fNtracks),
 262 fxOverX0Pipe(tracker.fxOverX0Pipe),
 263 fxTimesRhoPipe(tracker.fxTimesRhoPipe),
 264 fxOverX0PipeTrks(0),
 265 fxTimesRhoPipeTrks(0),
 266 fxOverX0ShieldTrks(0),
 267 fxTimesRhoShieldTrks(0),
 268 fxOverX0LayerTrks(0),
 269 fxTimesRhoLayerTrks(0),
 270 fDebugStreamer(tracker.fDebugStreamer),
 271 fITSChannelStatus(tracker.fITSChannelStatus),
 272 fkDetTypeRec(tracker.fkDetTypeRec),
 273 fPlaneEff(tracker.fPlaneEff) {
 274   //Copy constructor
 275   Int_t i;
 276   for(i=0;i<4;i++) {
 277     fSPDdetzcentre[i]=tracker.fSPDdetzcentre[i];
 278   }
 279   for(i=0;i<6;i++) {
 280     fxOverX0Layer[i]=tracker.fxOverX0Layer[i];
 281     fxTimesRhoLayer[i]=tracker.fxTimesRhoLayer[i];
 282   }
 283   for(i=0;i<2;i++) {
 284     fxOverX0Shield[i]=tracker.fxOverX0Shield[i];
 285     fxTimesRhoShield[i]=tracker.fxTimesRhoShield[i];
 286   }
 287 }
 288 //------------------------------------------------------------------------
 289 AliITStrackerMI & AliITStrackerMI::operator=(const AliITStrackerMI &tracker){
 290   //Assignment operator
 291   this->~AliITStrackerMI();
 292   new(this) AliITStrackerMI(tracker);
 293   return *this;
 294 }
 295 //------------------------------------------------------------------------
 296 AliITStrackerMI::~AliITStrackerMI()
 297 {
 298   //
 299   //destructor
 300   //
 301   if (fCoefficients) delete [] fCoefficients;
 302   DeleteTrksMaterialLUT();
 303   if (fDebugStreamer) {
 304     //fDebugStreamer->Close();
 305     delete fDebugStreamer;
 306   }
 307   if(fITSChannelStatus) delete fITSChannelStatus;
 308   if(fPlaneEff) delete fPlaneEff;
 309 }
 310 //------------------------------------------------------------------------
 311 void AliITStrackerMI::ReadBadFromDetTypeRec() {
 312   //--------------------------------------------------------------------
 313   //This function read ITS bad detectors, chips, channels from AliITSDetTypeRec
 314   //i.e. from OCDB
 315   //--------------------------------------------------------------------
 316
 317   if(!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return;
 318
 319   Info("ReadBadFromDetTypeRec","Reading info about bad ITS detectors and channels");
 320
 321   if(!fkDetTypeRec) Error("ReadBadFromDetTypeRec","AliITSDetTypeRec nof found!\n");
 322
 323   // ITS channels map
 324   if(fITSChannelStatus) delete fITSChannelStatus;
 325   fITSChannelStatus = new AliITSChannelStatus(fkDetTypeRec);
 326
 327   // ITS detectors and chips
 328   Int_t i=0,j=0,k=0,ndet=0;
 329   for (i=1; i<AliITSgeomTGeo::GetNLayers()+1; i++) {
 330     Int_t nBadDetsPerLayer=0;
 331     ndet=AliITSgeomTGeo::GetNDetectors(i);
 332     for (j=1; j<AliITSgeomTGeo::GetNLadders(i)+1; j++) {
 333       for (k=1; k<ndet+1; k++) {
 334         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
 335         det.ReadBadDetectorAndChips(i-1,(j-1)*ndet + k-1,fkDetTypeRec);
 336         if(det.IsBad()) {nBadDetsPerLayer++;}
 337       } // end loop on detectors
 338     } // end loop on ladders
 339     Info("ReadBadFromDetTypeRec",Form("Layer %d: %d bad out of %d",i-1,nBadDetsPerLayer,ndet*AliITSgeomTGeo::GetNLadders(i)));
 340   } // end loop on layers
 341
 342   return;
 343 }
 344 //------------------------------------------------------------------------
 345 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
 346   //--------------------------------------------------------------------
 347   //This function loads ITS clusters
 348   //--------------------------------------------------------------------
 349   TBranch *branch=cTree->GetBranch("ITSRecPoints");
 350   if (!branch) {
 351     Error("LoadClusters"," can't get the branch !\n");
 352     return 1;
 353   }
 354
 355   static TClonesArray dummy("AliITSRecPoint",10000), *clusters=&dummy;
 356   branch->SetAddress(&clusters);
 357
 358   Int_t i=0,j=0,ndet=0;
 359   Int_t detector=0;
 360   for (i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 361     ndet=fgLayers[i].GetNdetectors();
 362     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 363     for (; j<jmax; j++) {
 364       if (!cTree->GetEvent(j)) continue;
 365       Int_t ncl=clusters->GetEntriesFast();
 366       SignDeltas(clusters,GetZ());
 367
 368       while (ncl--) {
 369         AliITSRecPoint *c=(AliITSRecPoint*)clusters->UncheckedAt(ncl);
 370         detector=c->GetDetectorIndex();
 371
 372         if (!c->Misalign()) AliWarning("Can't misalign this cluster !");
 373
 374         fgLayers[i].InsertCluster(new AliITSRecPoint(*c));
 375       }
 376       clusters->Clear();
 377       // add dead zone "virtual" cluster in SPD, if there is a cluster within
 378       // zwindow cm from the dead zone
 379       if (i<2 && AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
 380         for (Float_t xdead = 0; xdead < AliITSRecoParam::GetSPDdetxlength(); xdead += (i+1.)*AliITSReconstructor::GetRecoParam()->GetXPassDeadZoneHits()) {
 381           Int_t lab[4]   = {0,0,0,detector};
 382           Int_t info[3]  = {0,0,i};
 383           Float_t q      = 0.; // this identifies virtual clusters
 384           Float_t hit[5] = {xdead,
 385                             0.,
 386                             AliITSReconstructor::GetRecoParam()->GetSigmaXDeadZoneHit2(),
 387                             AliITSReconstructor::GetRecoParam()->GetSigmaZDeadZoneHit2(),
 388                             q};
 389           Bool_t local   = kTRUE;
 390           Double_t zwindow = AliITSReconstructor::GetRecoParam()->GetZWindowDeadZone();
 391           hit[1] = fSPDdetzcentre[0]+0.5*AliITSRecoParam::GetSPDdetzlength();
 392           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 393             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 394           hit[1] = fSPDdetzcentre[1]-0.5*AliITSRecoParam::GetSPDdetzlength();
 395           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 396             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 397           hit[1] = fSPDdetzcentre[1]+0.5*AliITSRecoParam::GetSPDdetzlength();
 398           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 399             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 400           hit[1] = fSPDdetzcentre[2]-0.5*AliITSRecoParam::GetSPDdetzlength();
 401           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 402             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 403           hit[1] = fSPDdetzcentre[2]+0.5*AliITSRecoParam::GetSPDdetzlength();
 404           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 405             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 406           hit[1] = fSPDdetzcentre[3]-0.5*AliITSRecoParam::GetSPDdetzlength();
 407           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<zwindow)
 408             fgLayers[i].InsertCluster(new AliITSRecPoint(lab,hit,info,local));
 409         }
 410       } // "virtual" clusters in SPD
 411
 412     }
 413     //
 414     fgLayers[i].ResetRoad(); //road defined by the cluster density
 415     fgLayers[i].SortClusters();
 416   }
 417
 418   // check whether we have to skip some layers
 419   SetForceSkippingOfLayer();
 420
 421   return 0;
 422 }
 423 //------------------------------------------------------------------------
 424 void AliITStrackerMI::UnloadClusters() {
 425   //--------------------------------------------------------------------
 426   //This function unloads ITS clusters
 427   //--------------------------------------------------------------------
 428   for (Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) fgLayers[i].ResetClusters();
 429 }
 430 //------------------------------------------------------------------------
 431 void AliITStrackerMI::FillClusterArray(TObjArray* array) const {
 432   //--------------------------------------------------------------------
 433   // Publishes all pointers to clusters known to the tracker into the
 434   // passed object array.
 435   // The ownership is not transfered - the caller is not expected to delete
 436   // the clusters.
 437   //--------------------------------------------------------------------
 438
 439   for(Int_t i=0; i<AliITSgeomTGeo::GetNLayers(); i++) {
 440     for(Int_t icl=0; icl<fgLayers[i].GetNumberOfClusters(); icl++) {
 441       AliCluster *cl = (AliCluster*)fgLayers[i].GetCluster(icl);
 442       array->AddLast(cl);
 443     }
 444   }
 445
 446   return;
 447 }
 448 //------------------------------------------------------------------------
 449 Int_t AliITStrackerMI::CorrectForTPCtoITSDeadZoneMaterial(AliITStrackMI *t) {
 450   //--------------------------------------------------------------------
 451   // Correction for the material between the TPC and the ITS
 452   //--------------------------------------------------------------------
 453   if (t->GetX() > AliITSRecoParam::Getriw()) {   // inward direction
 454       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw(),1)) return 0;// TPC inner wall
 455       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1)) return 0;// TPC central drum
 456       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))  return 0;// ITS screen
 457   } else if (t->GetX() < AliITSRecoParam::Getrs()) {  // outward direction
 458       if (!t->PropagateToTGeo(AliITSRecoParam::Getrs(),1))        return 0;// ITS screen
 459       if (!t->PropagateToTGeo(AliITSRecoParam::Getrcd(),1))       return 0;// TPC central drum
 460       if (!t->PropagateToTGeo(AliITSRecoParam::Getriw()+0.001,1)) return 0;// TPC inner wall
 461   } else {
 462     printf("CorrectForTPCtoITSDeadZoneMaterial: Track is already in the dead zone !\n");
 463     return 0;
 464   }
 465
 466   return 1;
 467 }
 468 //------------------------------------------------------------------------
 469 Int_t AliITStrackerMI::Clusters2Tracks(AliESDEvent *event) {
 470   //--------------------------------------------------------------------
 471   // This functions reconstructs ITS tracks
 472   // The clusters must be already loaded !
 473   //--------------------------------------------------------------------
 474
 475   AliDebug(2,Form("SKIPPING %d %d %d %d %d %d",ForceSkippingOfLayer(0),ForceSkippingOfLayer(1),ForceSkippingOfLayer(2),ForceSkippingOfLayer(3),ForceSkippingOfLayer(4),ForceSkippingOfLayer(5)));
 476
 477   fTrackingPhase="Clusters2Tracks";
 478
 479   TObjArray itsTracks(15000);
 480   fOriginal.Clear();
 481   fEsd = event;         // store pointer to the esd
 482
 483   // temporary (for cosmics)
 484   if(event->GetVertex()) {
 485     TString title = event->GetVertex()->GetTitle();
 486     if(title.Contains("cosmics")) {
 487       Double_t xyz[3]={GetX(),GetY(),GetZ()};
 488       Double_t exyz[3]={0.1,0.1,0.1};
 489       SetVertex(xyz,exyz);
 490     }
 491   }
 492   // temporary
 493   Int_t noesd = 0;
 494   {/* Read ESD tracks */
 495     Double_t pimass = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 496     Int_t nentr=event->GetNumberOfTracks();
 497     noesd=nentr;
 498     //    Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 499     while (nentr--) {
 500       AliESDtrack *esd=event->GetTrack(nentr);
 501       //  ---- for debugging:
 502       //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); }
 503
 504       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 505       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 506       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 507       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
 508       AliITStrackMI *t=0;
 509       try {
 510         t=new AliITStrackMI(*esd);
 511       } catch (const Char_t *msg) {
 512         //Warning("Clusters2Tracks",msg);
 513         delete t;
 514         continue;
 515       }
 516       t->GetDZ(GetX(),GetY(),GetZ(),t->GetDP());              //I.B.
 517       Double_t vdist = TMath::Sqrt(t->GetD(0)*t->GetD(0)+t->GetD(1)*t->GetD(1));
 518
 519
 520       // look at the ESD mass hypothesys !
 521       if (t->GetMass()<0.9*pimass) t->SetMass(pimass);
 522       // write expected q
 523       t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 524
 525       if (esd->GetV0Index(0)>0 && t->GetD(0)<AliITSReconstructor::GetRecoParam()->GetMaxDforV0dghtrForProlongation()){
 526         //track - can be  V0 according to TPC
 527       } else {
 528         if (TMath::Abs(t->GetD(0))>AliITSReconstructor::GetRecoParam()->GetMaxDForProlongation()) {
 529           delete t;
 530           continue;
 531         }
 532         if (TMath::Abs(vdist)>AliITSReconstructor::GetRecoParam()->GetMaxDZForProlongation()) {
 533           delete t;
 534           continue;
 535         }
 536         if (t->Pt()<AliITSReconstructor::GetRecoParam()->GetMinPtForProlongation()) {
 537           delete t;
 538           continue;
 539         }
 540         if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 541           delete t;
 542           continue;
 543         }
 544       }
 545       t->SetReconstructed(kFALSE);
 546       itsTracks.AddLast(t);
 547       fOriginal.AddLast(t);
 548     }
 549   } /* End Read ESD tracks */
 550
 551   itsTracks.Sort();
 552   fOriginal.Sort();
 553   Int_t nentr=itsTracks.GetEntriesFast();
 554   fTrackHypothesys.Expand(nentr);
 555   fBestHypothesys.Expand(nentr);
 556   MakeCoefficients(nentr);
 557   if(fUseTGeo==3 || fUseTGeo==4) MakeTrksMaterialLUT(event->GetNumberOfTracks());
 558   Int_t ntrk=0;
 559   // THE TWO TRACKING PASSES
 560   for (fPass=0; fPass<2; fPass++) {
 561      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 562      for (fCurrentEsdTrack=0; fCurrentEsdTrack<nentr; fCurrentEsdTrack++) {
 563        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(fCurrentEsdTrack);
 564        if (t==0) continue;              //this track has been already tracked
 565        //cout<<"========== "<<fPass<<"    "<<fCurrentEsdTrack<<" =========\n";
 566        if (t->GetReconstructed()&&(t->GetNUsed()<1.5)) continue;  //this track was  already  "succesfully" reconstructed
 567        Float_t dz[2]; t->GetDZ(GetX(),GetY(),GetZ(),dz);              //I.B.
 568        if (fConstraint[fPass]) {
 569          if (TMath::Abs(dz[0])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint() ||
 570              TMath::Abs(dz[1])>AliITSReconstructor::GetRecoParam()->GetMaxDZToUseConstraint()) continue;
 571        }
 572
 573        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 574        AliDebug(2,Form("LABEL %d pass %d",tpcLabel,fPass));
 575        fI = 6;
 576        ResetTrackToFollow(*t);
 577        ResetBestTrack();
 578
 579        FollowProlongationTree(t,fCurrentEsdTrack,fConstraint[fPass]);
 580
 581
 582        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
 583        //
 584        AliITStrackMI *besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
 585        if (!besttrack) continue;
 586        besttrack->SetLabel(tpcLabel);
 587        //       besttrack->CookdEdx();
 588        CookdEdx(besttrack);
 589        besttrack->SetFakeRatio(1.);
 590        CookLabel(besttrack,0.); //For comparison only
 591        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
 592
 593        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain
 594
 595        t->SetReconstructed(kTRUE);
 596        ntrk++;
 597        AliDebug(2,Form("TRACK! (label %d) ncls %d",besttrack->GetLabel(),besttrack->GetNumberOfClusters()));
 598      }
 599      GetBestHypothesysMIP(itsTracks);
 600   } // end loop on the two tracking passes
 601
 602   if(event->GetNumberOfV0s()>0) AliITSV0Finder::UpdateTPCV0(event,this);
 603   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::FindV02(event,this);
 604   fAfterV0 = kTRUE;
 605   //
 606   itsTracks.Delete();
 607   //
 608   Int_t entries = fTrackHypothesys.GetEntriesFast();
 609   for (Int_t ientry=0; ientry<entries; ientry++) {
 610     TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
 611     if (array) array->Delete();
 612     delete fTrackHypothesys.RemoveAt(ientry);
 613   }
 614
 615   fTrackHypothesys.Delete();
 616   fBestHypothesys.Delete();
 617   fOriginal.Clear();
 618   delete [] fCoefficients;
 619   fCoefficients=0;
 620   DeleteTrksMaterialLUT();
 621
 622   AliInfo(Form("Number of prolonged tracks: %d out of %d ESD tracks",ntrk,noesd));
 623
 624   fTrackingPhase="Default";
 625
 626   return 0;
 627 }
 628 //------------------------------------------------------------------------
 629 Int_t AliITStrackerMI::PropagateBack(AliESDEvent *event) {
 630   //--------------------------------------------------------------------
 631   // This functions propagates reconstructed ITS tracks back
 632   // The clusters must be loaded !
 633   //--------------------------------------------------------------------
 634   fTrackingPhase="PropagateBack";
 635   Int_t nentr=event->GetNumberOfTracks();
 636   //  Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 637
 638   Int_t ntrk=0;
 639   for (Int_t i=0; i<nentr; i++) {
 640      AliESDtrack *esd=event->GetTrack(i);
 641
 642      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 643      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 644
 645      AliITStrackMI *t=0;
 646      try {
 647         t=new AliITStrackMI(*esd);
 648      } catch (const Char_t *msg) {
 649        //Warning("PropagateBack",msg);
 650         delete t;
 651         continue;
 652      }
 653      t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 654
 655      ResetTrackToFollow(*t);
 656
 657      /*
 658      // propagate to vertex [SR, GSI 17.02.2003]
 659      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
 660      if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 661        if (fTrackToFollow.PropagateToVertex(event->GetVertex()))
 662          fTrackToFollow.StartTimeIntegral();
 663        // from vertex to outside pipe
 664        CorrectForPipeMaterial(&fTrackToFollow,"outward");
 665        }*/
 666      // Start time integral and add distance from current position to vertex
 667      Double_t xyzTrk[3],xyzVtx[3]={GetX(),GetY(),GetZ()};
 668      fTrackToFollow.GetXYZ(xyzTrk);
 669      Double_t dst2 = 0.;
 670      for (Int_t icoord=0; icoord<3; icoord++) {
 671        Double_t di = xyzTrk[icoord] - xyzVtx[icoord];
 672        dst2 += di*di;
 673      }
 674      fTrackToFollow.StartTimeIntegral();
 675      fTrackToFollow.AddTimeStep(TMath::Sqrt(dst2));
 676
 677      fTrackToFollow.ResetCovariance(10.); fTrackToFollow.ResetClusters();
 678      if (RefitAt(AliITSRecoParam::GetrInsideITSscreen(),&fTrackToFollow,t)) {
 679         if (!CorrectForTPCtoITSDeadZoneMaterial(&fTrackToFollow)) {
 680           delete t;
 681           continue;
 682         }
 683         fTrackToFollow.SetLabel(t->GetLabel());
 684         //fTrackToFollow.CookdEdx();
 685         CookLabel(&fTrackToFollow,0.); //For comparison only
 686         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 687         //UseClusters(&fTrackToFollow);
 688         ntrk++;
 689      }
 690      delete t;
 691   }
 692
 693   AliInfo(Form("Number of back propagated ITS tracks: %d out of %d ESD tracks",ntrk,nentr));
 694
 695   fTrackingPhase="Default";
 696
 697   return 0;
 698 }
 699 //------------------------------------------------------------------------
 700 Int_t AliITStrackerMI::RefitInward(AliESDEvent *event) {
 701   //--------------------------------------------------------------------
 702   // This functions refits ITS tracks using the
 703   // "inward propagated" TPC tracks
 704   // The clusters must be loaded !
 705   //--------------------------------------------------------------------
 706   fTrackingPhase="RefitInward";
 707
 708   if(AliITSReconstructor::GetRecoParam()->GetFindV0s()) AliITSV0Finder::RefitV02(event,this);
 709
 710   Int_t nentr=event->GetNumberOfTracks();
 711   //  Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 712
 713   Int_t ntrk=0;
 714   for (Int_t i=0; i<nentr; i++) {
 715     AliESDtrack *esd=event->GetTrack(i);
 716
 717     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 718     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 719     if (esd->GetStatus()&AliESDtrack::kTPCout)
 720       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 721
 722     AliITStrackMI *t=0;
 723     try {
 724         t=new AliITStrackMI(*esd);
 725     } catch (const Char_t *msg) {
 726       //Warning("RefitInward",msg);
 727         delete t;
 728         continue;
 729     }
 730     t->SetExpQ(TMath::Max(0.8*t->GetESDtrack()->GetTPCsignal(),30.));
 731     if (!CorrectForTPCtoITSDeadZoneMaterial(t)) {
 732        delete t;
 733        continue;
 734     }
 735
 736     ResetTrackToFollow(*t);
 737     fTrackToFollow.ResetClusters();
 738
 739     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
 740       fTrackToFollow.ResetCovariance(10.);
 741
 742     //Refitting...
 743     Bool_t pe=(AliITSReconstructor::GetRecoParam()->GetComputePlaneEff() &&
 744                AliITSReconstructor::GetRecoParam()->GetIPlanePlaneEff()>=0);
 745
 746     AliDebug(2,Form("Refit LABEL %d  %d",t->GetLabel(),t->GetNumberOfClusters()));
 747     if (RefitAt(AliITSRecoParam::GetrInsideSPD1(),&fTrackToFollow,t,kTRUE,pe)) {
 748        AliDebug(2,"  refit OK");
 749        fTrackToFollow.SetLabel(t->GetLabel());
 750        //       fTrackToFollow.CookdEdx();
 751        CookdEdx(&fTrackToFollow);
 752
 753        CookLabel(&fTrackToFollow,0.0); //For comparison only
 754
 755        //The beam pipe
 756        if (CorrectForPipeMaterial(&fTrackToFollow,"inward")) {
 757          fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
 758          AliESDtrack  *esdTrack =fTrackToFollow.GetESDtrack();
 759          //printf("                                       %d\n",esdTrack->GetITSModuleIndex(0));
 760          //esdTrack->UpdateTrackParams(&fTrackToFollow,AliESDtrack::kITSrefit); //original line
 761          Double_t r[3]={0.,0.,0.};
 762          Double_t maxD=3.;
 763          esdTrack->RelateToVertex(event->GetVertex(),GetBz(r),maxD);
 764          ntrk++;
 765        }
 766     }
 767     delete t;
 768   }
 769
 770   AliInfo(Form("Number of refitted tracks: %d out of %d ESD tracks",ntrk,nentr));
 771
 772   fTrackingPhase="Default";
 773
 774   return 0;
 775 }
 776 //------------------------------------------------------------------------
 777 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 778   //--------------------------------------------------------------------
 779   //       Return pointer to a given cluster
 780   //--------------------------------------------------------------------
 781   Int_t l=(index & 0xf0000000) >> 28;
 782   Int_t c=(index & 0x0fffffff) >> 00;
 783   return fgLayers[l].GetCluster(c);
 784 }
 785 //------------------------------------------------------------------------
 786 Bool_t AliITStrackerMI::GetTrackPoint(Int_t index, AliTrackPoint& p) const {
 787   //--------------------------------------------------------------------
 788   // Get track space point with index i
 789   //--------------------------------------------------------------------
 790
 791   Int_t l=(index & 0xf0000000) >> 28;
 792   Int_t c=(index & 0x0fffffff) >> 00;
 793   AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 794   Int_t idet = cl->GetDetectorIndex();
 795
 796   Float_t xyz[3];
 797   Float_t cov[6];
 798   cl->GetGlobalXYZ(xyz);
 799   cl->GetGlobalCov(cov);
 800   p.SetXYZ(xyz, cov);
 801   p.SetCharge(cl->GetQ());
 802   p.SetDriftTime(cl->GetDriftTime());
 803   p.SetChargeRatio(cl->GetChargeRatio());
 804   p.SetClusterType(cl->GetClusterType());
 805   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 806   switch (l) {
 807   case 0:
 808     iLayer = AliGeomManager::kSPD1;
 809     break;
 810   case 1:
 811     iLayer = AliGeomManager::kSPD2;
 812     break;
 813   case 2:
 814     iLayer = AliGeomManager::kSDD1;
 815     break;
 816   case 3:
 817     iLayer = AliGeomManager::kSDD2;
 818     break;
 819   case 4:
 820     iLayer = AliGeomManager::kSSD1;
 821     break;
 822   case 5:
 823     iLayer = AliGeomManager::kSSD2;
 824     break;
 825   default:
 826     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 827     break;
 828   };
 829   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 830   p.SetVolumeID((UShort_t)volid);
 831   return kTRUE;
 832 }
 833 //------------------------------------------------------------------------
 834 Bool_t AliITStrackerMI::GetTrackPointTrackingError(Int_t index,
 835                         AliTrackPoint& p, const AliESDtrack *t) {
 836   //--------------------------------------------------------------------
 837   // Get track space point with index i
 838   // (assign error estimated during the tracking)
 839   //--------------------------------------------------------------------
 840
 841   Int_t l=(index & 0xf0000000) >> 28;
 842   Int_t c=(index & 0x0fffffff) >> 00;
 843   const AliITSRecPoint *cl = fgLayers[l].GetCluster(c);
 844   Int_t idet = cl->GetDetectorIndex();
 845
 846   const AliITSdetector &det=fgLayers[l].GetDetector(idet);
 847
 848   // tgphi and tglambda of the track in tracking frame with alpha=det.GetPhi
 849   Float_t detxy[2];
 850   detxy[0] = det.GetR()*TMath::Cos(det.GetPhi());
 851   detxy[1] = det.GetR()*TMath::Sin(det.GetPhi());
 852   Double_t alpha = t->GetAlpha();
 853   Double_t xdetintrackframe = detxy[0]*TMath::Cos(alpha)+detxy[1]*TMath::Sin(alpha);
 854   Float_t phi = TMath::ASin(t->GetSnpAt(xdetintrackframe,GetBz()));
 855   phi += alpha-det.GetPhi();
 856   Float_t tgphi = TMath::Tan(phi);
 857
 858   Float_t tgl = t->GetTgl(); // tgl about const along track
 859   Float_t expQ = TMath::Max(0.8*t->GetTPCsignal(),30.);
 860
 861   Float_t errtrky,errtrkz,covyz;
 862   Bool_t addMisalErr=kFALSE;
 863   AliITSClusterParam::GetError(l,cl,tgl,tgphi,expQ,errtrky,errtrkz,covyz,addMisalErr);
 864
 865   Float_t xyz[3];
 866   Float_t cov[6];
 867   cl->GetGlobalXYZ(xyz);
 868   //  cl->GetGlobalCov(cov);
 869   Float_t pos[3] = {0.,0.,0.};
 870   AliCluster tmpcl((UShort_t)cl->GetVolumeId(),pos[0],pos[1],pos[2],errtrky*errtrky,errtrkz*errtrkz,covyz);
 871   tmpcl.GetGlobalCov(cov);
 872
 873   p.SetXYZ(xyz, cov);
 874   p.SetCharge(cl->GetQ());
 875   p.SetDriftTime(cl->GetDriftTime());
 876   p.SetChargeRatio(cl->GetChargeRatio());
 877   p.SetClusterType(cl->GetClusterType());
 878
 879   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
 880   switch (l) {
 881   case 0:
 882     iLayer = AliGeomManager::kSPD1;
 883     break;
 884   case 1:
 885     iLayer = AliGeomManager::kSPD2;
 886     break;
 887   case 2:
 888     iLayer = AliGeomManager::kSDD1;
 889     break;
 890   case 3:
 891     iLayer = AliGeomManager::kSDD2;
 892     break;
 893   case 4:
 894     iLayer = AliGeomManager::kSSD1;
 895     break;
 896   case 5:
 897     iLayer = AliGeomManager::kSSD2;
 898     break;
 899   default:
 900     AliWarning(Form("Wrong layer index in ITS (%d) !",l));
 901     break;
 902   };
 903   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,idet);
 904
 905   p.SetVolumeID((UShort_t)volid);
 906   return kTRUE;
 907 }
 908 //------------------------------------------------------------------------
 909 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain)
 910 {
 911   //--------------------------------------------------------------------
 912   // Follow prolongation tree
 913   //--------------------------------------------------------------------
 914   //
 915   Double_t xyzVtx[]={GetX(),GetY(),GetZ()};
 916   Double_t ersVtx[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 917
 918
 919   AliESDtrack * esd = otrack->GetESDtrack();
 920   if (esd->GetV0Index(0)>0) {
 921     // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
 922     //                      mapping of ESD track is different as ITS track in Containers
 923     //                      Need something more stable
 924     //                      Indexes are set back again to the ESD track indexes in UpdateTPCV0
 925     for (Int_t i=0;i<3;i++){
 926       Int_t  index = esd->GetV0Index(i);
 927       if (index==0) break;
 928       AliESDv0 * vertex = fEsd->GetV0(index);
 929       if (vertex->GetStatus()<0) continue;     // rejected V0
 930       //
 931       if (esd->GetSign()>0) {
 932         vertex->SetIndex(0,esdindex);
 933       } else {
 934         vertex->SetIndex(1,esdindex);
 935       }
 936     }
 937   }
 938   TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
 939   if (!bestarray){
 940     bestarray = new TObjArray(5);
 941     fBestHypothesys.AddAt(bestarray,esdindex);
 942   }
 943
 944   //
 945   //setup tree of the prolongations
 946   //
 947   static AliITStrackMI tracks[7][100];
 948   AliITStrackMI *currenttrack;
 949   static AliITStrackMI currenttrack1;
 950   static AliITStrackMI currenttrack2;
 951   static AliITStrackMI backuptrack;
 952   Int_t ntracks[7];
 953   Int_t nindexes[7][100];
 954   Float_t normalizedchi2[100];
 955   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
 956   otrack->SetNSkipped(0);
 957   new (&(tracks[6][0])) AliITStrackMI(*otrack);
 958   ntracks[6]=1;
 959   for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
 960   Int_t modstatus = 1; // found
 961   Float_t xloc,zloc;
 962   //
 963   //
 964   // follow prolongations
 965   for (Int_t ilayer=5; ilayer>=0; ilayer--) {
 966     AliDebug(2,Form("FollowProlongationTree: layer %d",ilayer));
 967     fI = ilayer;
 968     //
 969     AliITSlayer &layer=fgLayers[ilayer];
 970     Double_t r = layer.GetR();
 971     ntracks[ilayer]=0;
 972     //
 973     //
 974     Int_t nskipped=0;
 975     Float_t nused =0;
 976     for (Int_t itrack =0; itrack<ntracks[ilayer+1]; itrack++) {
 977       //set current track
 978       if (ntracks[ilayer]>=100) break;
 979       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0) nskipped++;
 980       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2.) nused++;
 981       if (ntracks[ilayer]>15+ilayer){
 982         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNSkipped()>0 && nskipped>4+ilayer) continue;
 983         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].GetNUsed()>2. && nused>3) continue;
 984       }
 985
 986       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
 987
 988       // material between SSD and SDD, SDD and SPD
 989       if (ilayer==3)
 990         if(!CorrectForShieldMaterial(&currenttrack1,"SDD","inward")) continue;
 991       if (ilayer==1)
 992         if(!CorrectForShieldMaterial(&currenttrack1,"SPD","inward")) continue;
 993
 994       // detector number
 995       Double_t phi,z;
 996       if (!currenttrack1.GetPhiZat(r,phi,z)) continue;
 997       Int_t idet=layer.FindDetectorIndex(phi,z);
 998
 999       Double_t trackGlobXYZ1[3];
1000       if (!currenttrack1.GetXYZ(trackGlobXYZ1)) continue;
1001
1002       // Get the budget to the primary vertex for the current track being prolonged
1003       Double_t budgetToPrimVertex = GetEffectiveThickness();
1004
1005       // check if we allow a prolongation without point
1006       Int_t skip = CheckSkipLayer(&currenttrack1,ilayer,idet);
1007       if (skip) {
1008         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1009         // propagate to the layer radius
1010         Double_t xToGo; if (!vtrack->GetLocalXat(r,xToGo)) continue;
1011         if(!vtrack->Propagate(xToGo)) continue;
1012         // apply correction for material of the current layer
1013         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1014         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1015         vtrack->SetDeadZoneProbability(ilayer,1.); // no penalty for missing cluster
1016         vtrack->SetClIndex(ilayer,-1);
1017         modstatus = (skip==1 ? 3 : 4); // skipped : out in z
1018         if(LocalModuleCoord(ilayer,idet,vtrack,xloc,zloc)) { // local module coords
1019           vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1020         }
1021         if(constrain) vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1022         ntracks[ilayer]++;
1023         continue;
1024       }
1025
1026       // track outside layer acceptance in z
1027       if (idet<0) continue;
1028
1029       //propagate to the intersection with the detector plane
1030       const AliITSdetector &det=layer.GetDetector(idet);
1031       new(&currenttrack2)  AliITStrackMI(currenttrack1);
1032       if (!currenttrack1.Propagate(det.GetPhi(),det.GetR())) continue;
1033       if (!currenttrack2.Propagate(det.GetPhi(),det.GetR())) continue;
1034       currenttrack1.SetDetectorIndex(idet);
1035       currenttrack2.SetDetectorIndex(idet);
1036       if(!LocalModuleCoord(ilayer,idet,&currenttrack1,xloc,zloc)) continue; // local module coords
1037
1038       //***************
1039       // DEFINITION OF SEARCH ROAD AND CLUSTERS SELECTION
1040       //
1041       // road in global (rphi,z) [i.e. in tracking ref. system]
1042       Double_t zmin,zmax,ymin,ymax;
1043       if (!ComputeRoad(&currenttrack1,ilayer,idet,zmin,zmax,ymin,ymax)) continue;
1044
1045       // select clusters in road
1046       layer.SelectClusters(zmin,zmax,ymin,ymax);
1047       //********************
1048
1049       // Define criteria for track-cluster association
1050       Double_t msz = currenttrack1.GetSigmaZ2() +
1051         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1052         AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1053         AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
1054       Double_t msy = currenttrack1.GetSigmaY2() +
1055         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1056         AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1057         AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
1058       if (constrain) {
1059         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
1060         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
1061       }  else {
1062         msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
1063         msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
1064       }
1065       msz = 1./msz; // 1/RoadZ^2
1066       msy = 1./msy; // 1/RoadY^2
1067
1068       //
1069       //
1070       // LOOP OVER ALL POSSIBLE TRACK PROLONGATIONS ON THIS LAYER
1071       //
1072       const AliITSRecPoint *cl=0;
1073       Int_t clidx=-1;
1074       Double_t chi2trkcl=AliITSReconstructor::GetRecoParam()->GetMaxChi2(); // init with big value
1075       Bool_t deadzoneSPD=kFALSE;
1076       currenttrack = &currenttrack1;
1077
1078       // check if the road contains a dead zone
1079       Bool_t noClusters = kFALSE;
1080       if (!layer.GetNextCluster(clidx,kTRUE)) noClusters=kTRUE;
1081       if (noClusters) AliDebug(2,"no clusters in road");
1082       Double_t dz=0.5*(zmax-zmin);
1083       Double_t dy=0.5*(ymax-ymin);
1084       Int_t dead = CheckDeadZone(&currenttrack1,ilayer,idet,dz,dy,noClusters);
1085       if(dead) AliDebug(2,Form("DEAD (%d)\n",dead));
1086       // create a prolongation without clusters (check also if there are no clusters in the road)
1087       if (dead ||
1088           (noClusters &&
1089            AliITSReconstructor::GetRecoParam()->GetAllowProlongationWithEmptyRoad())) {
1090         AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1091         updatetrack->SetClIndex(ilayer,-1);
1092         if (dead==0) {
1093           modstatus = 5; // no cls in road
1094         } else if (dead==1) {
1095           modstatus = 7; // holes in z in SPD
1096         } else if (dead==2 || dead==3 || dead==4) {
1097           modstatus = 2; // dead from OCDB
1098         }
1099         updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1100         // apply correction for material of the current layer
1101         CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1102         if (constrain) { // apply vertex constrain
1103           updatetrack->SetConstrain(constrain);
1104           Bool_t isPrim = kTRUE;
1105           if (ilayer<4) { // check that it's close to the vertex
1106             updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1107             if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1108                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1109                 TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1110                 AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1111           }
1112           if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1113         }
1114         updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1115         if (dead) {
1116           if (dead==1) { // dead zone at z=0,+-7cm in SPD
1117             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1118             deadzoneSPD=kTRUE;
1119           } else if (dead==2 || dead==3) { // dead module or chip from OCDB
1120             updatetrack->SetDeadZoneProbability(ilayer,1.);
1121           } else if (dead==4) { // at least a single dead channel from OCDB
1122             updatetrack->SetDeadZoneProbability(ilayer,0.);
1123           }
1124         }
1125         ntracks[ilayer]++;
1126       }
1127
1128       clidx=-1;
1129       // loop over clusters in the road
1130       while ((cl=layer.GetNextCluster(clidx))!=0) {
1131         if (ntracks[ilayer]>95) break; //space for skipped clusters
1132         Bool_t changedet =kFALSE;
1133         if (cl->GetQ()==0 && deadzoneSPD==kTRUE) continue;
1134         Int_t idetc=cl->GetDetectorIndex();
1135
1136         if (currenttrack->GetDetectorIndex()==idetc) { // track already on the cluster's detector
1137           // take into account misalignment (bring track to real detector plane)
1138           Double_t xTrOrig = currenttrack->GetX();
1139           if (!currenttrack->Propagate(xTrOrig+cl->GetX())) continue;
1140           // a first cut on track-cluster distance
1141           if ( (currenttrack->GetZ()-cl->GetZ())*(currenttrack->GetZ()-cl->GetZ())*msz +
1142                (currenttrack->GetY()-cl->GetY())*(currenttrack->GetY()-cl->GetY())*msy > 1. )
1143             {  // cluster not associated to track
1144               AliDebug(2,"not associated");
1145               continue;
1146             }
1147           // bring track back to ideal detector plane
1148           if (!currenttrack->Propagate(xTrOrig)) continue;
1149         } else {                                      // have to move track to cluster's detector
1150           const AliITSdetector &detc=layer.GetDetector(idetc);
1151           // a first cut on track-cluster distance
1152           Double_t y;
1153           if (!currenttrack2.GetProlongationFast(detc.GetPhi(),detc.GetR()+cl->GetX(),y,z)) continue;
1154           if ( (z-cl->GetZ())*(z-cl->GetZ())*msz +
1155                (y-cl->GetY())*(y-cl->GetY())*msy > 1. )
1156             continue; // cluster not associated to track
1157           //
1158           new (&backuptrack) AliITStrackMI(currenttrack2);
1159           changedet = kTRUE;
1160           currenttrack =&currenttrack2;
1161           if (!currenttrack->Propagate(detc.GetPhi(),detc.GetR())) {
1162             new (currenttrack) AliITStrackMI(backuptrack);
1163             changedet = kFALSE;
1164             continue;
1165           }
1166           currenttrack->SetDetectorIndex(idetc);
1167           // Get again the budget to the primary vertex
1168           // for the current track being prolonged, if had to change detector
1169           //budgetToPrimVertex = GetEffectiveThickness();// not needed at the moment because anyway we take a mean material for this correction
1170         }
1171
1172         // calculate track-clusters chi2
1173         chi2trkcl = GetPredictedChi2MI(currenttrack,cl,ilayer);
1174         // chi2 cut
1175         AliDebug(2,Form("chi2 %f max %f",chi2trkcl,AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)));
1176         if (chi2trkcl < AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) {
1177           if (cl->GetQ()==0) deadzoneSPD=kTRUE; // only 1 prolongation with virtual cluster
1178           if (ntracks[ilayer]>=100) continue;
1179           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
1180           updatetrack->SetClIndex(ilayer,-1);
1181           if (changedet) new (&currenttrack2) AliITStrackMI(backuptrack);
1182
1183           if (cl->GetQ()!=0) { // real cluster
1184             if (!UpdateMI(updatetrack,cl,chi2trkcl,(ilayer<<28)+clidx)) {
1185               AliDebug(2,"update failed");
1186               continue;
1187             }
1188             updatetrack->SetSampledEdx(cl->GetQ(),ilayer-2);
1189             modstatus = 1; // found
1190           } else {             // virtual cluster in dead zone
1191             updatetrack->SetNDeadZone(updatetrack->GetNDeadZone()+1);
1192             updatetrack->SetDeadZoneProbability(ilayer,GetSPDDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2())));
1193             modstatus = 7; // holes in z in SPD
1194           }
1195
1196           if (changedet) {
1197             Float_t xlocnewdet,zlocnewdet;
1198             if(LocalModuleCoord(ilayer,idet,updatetrack,xlocnewdet,zlocnewdet)) { // local module coords
1199               updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xlocnewdet,zlocnewdet);
1200             }
1201           } else {
1202             updatetrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1203           }
1204           if (cl->IsUsed()) updatetrack->IncrementNUsed();
1205
1206           // apply correction for material of the current layer
1207           CorrectForLayerMaterial(updatetrack,ilayer,trackGlobXYZ1,"inward");
1208
1209           if (constrain) { // apply vertex constrain
1210             updatetrack->SetConstrain(constrain);
1211             Bool_t isPrim = kTRUE;
1212             if (ilayer<4) { // check that it's close to the vertex
1213               updatetrack->GetDZ(GetX(),GetY(),GetZ(),updatetrack->GetDP()); //I.B.
1214               if (TMath::Abs(updatetrack->GetD(0)/(1.+ilayer)) > // y
1215                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk() ||
1216                   TMath::Abs(updatetrack->GetD(1)/(1.+ilayer)) > // z
1217                   AliITSReconstructor::GetRecoParam()->GetMaxDZforPrimTrk()) isPrim=kFALSE;
1218             }
1219             if (isPrim) updatetrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1220           } //apply vertex constrain
1221           ntracks[ilayer]++;
1222         }  // create new hypothesis
1223         else {
1224           AliDebug(2,"chi2 too large");
1225         }
1226
1227       } // loop over possible prolongations
1228
1229       // allow one prolongation without clusters
1230       if (constrain && itrack<=1 && currenttrack1.GetNSkipped()==0 && deadzoneSPD==kFALSE && ntracks[ilayer]<100) {
1231         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1232         // apply correction for material of the current layer
1233         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1234         vtrack->SetClIndex(ilayer,-1);
1235         modstatus = 3; // skipped
1236         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1237         vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1238         vtrack->IncrementNSkipped();
1239         ntracks[ilayer]++;
1240       }
1241       /*
1242       // allow one prolongation without clusters for tracks with |tgl|>1.1
1243       if (constrain && itrack==0 && TMath::Abs(currenttrack1.GetTgl())>1.1) {  //big theta - for low flux
1244         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
1245         // apply correction for material of the current layer
1246         CorrectForLayerMaterial(vtrack,ilayer,trackGlobXYZ1,"inward");
1247         vtrack->SetClIndex(ilayer,-1);
1248         modstatus = 3; // skipped
1249         vtrack->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
1250         vtrack->Improve(budgetToPrimVertex,xyzVtx,ersVtx);
1251         vtrack->SetNDeadZone(vtrack->GetNDeadZone()+1);
1252         ntracks[ilayer]++;
1253       }
1254       */
1255
1256     } // loop over tracks in layer ilayer+1
1257
1258     //loop over track candidates for the current layer
1259     //
1260     //
1261     Int_t accepted=0;
1262
1263     Int_t golden=0;
1264     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
1265       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
1266       if (normalizedchi2[itrack] <
1267           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2ForGolden(ilayer)) golden++;
1268       if (ilayer>4) {
1269         accepted++;
1270       } else {
1271         if (constrain) { // constrain
1272           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2C(ilayer)+1)
1273             accepted++;
1274         } else {        // !constrain
1275           if (normalizedchi2[itrack]<AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonC(ilayer)+1)
1276             accepted++;
1277         }
1278       }
1279     }
1280     // sort tracks by increasing normalized chi2
1281     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
1282     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
1283     if (ntracks[ilayer]<golden+2+ilayer) ntracks[ilayer]=TMath::Min(golden+2+ilayer,accepted);
1284     if (ntracks[ilayer]>90) ntracks[ilayer]=90;
1285   } // end loop over layers
1286
1287
1288   //
1289   // Now select tracks to be kept
1290   //
1291   Int_t max = constrain ? 20 : 5;
1292
1293   // tracks that reach layer 0 (SPD inner)
1294   for (Int_t i=0; i<TMath::Min(max,ntracks[0]); i++) {
1295     AliITStrackMI & track= tracks[0][nindexes[0][i]];
1296     if (track.GetNumberOfClusters()<2) continue;
1297     if (!constrain && track.GetNormChi2(0) >
1298         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) {
1299       continue;
1300     }
1301     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1302   }
1303
1304   // tracks that reach layer 1 (SPD outer)
1305   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
1306     AliITStrackMI & track= tracks[1][nindexes[1][i]];
1307     if (track.GetNumberOfClusters()<4) continue;
1308     if (!constrain && track.GetNormChi2(1) >
1309         AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1310     if (constrain) track.IncrementNSkipped();
1311     if (!constrain) {
1312       track.SetD(0,track.GetD(GetX(),GetY()));
1313       track.SetNSkipped(track.GetNSkipped()+4./(4.+8.*TMath::Abs(track.GetD(0))));
1314       if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1315         track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1316       }
1317     }
1318     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1319   }
1320
1321   // tracks that reach layer 2 (SDD inner), only during non-constrained pass
1322   if (!constrain){
1323     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
1324       AliITStrackMI & track= tracks[2][nindexes[2][i]];
1325       if (track.GetNumberOfClusters()<3) continue;
1326       if (!constrain && track.GetNormChi2(2) >
1327           AliITSReconstructor::GetRecoParam()->GetMaxNormChi2NonCForHypothesis()) continue;
1328       if (constrain) track.SetNSkipped(track.GetNSkipped()+2);
1329       if (!constrain){
1330         track.SetD(0,track.GetD(GetX(),GetY()));
1331         track.SetNSkipped(track.GetNSkipped()+7./(7.+8.*TMath::Abs(track.GetD(0))));
1332         if (track.GetNumberOfClusters()+track.GetNDeadZone()+track.GetNSkipped()>6) {
1333           track.SetNSkipped(6-track.GetNumberOfClusters()+track.GetNDeadZone());
1334         }
1335       }
1336       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1337     }
1338   }
1339
1340   if (!constrain) {
1341     //
1342     // register best track of each layer - important for V0 finder
1343     //
1344     for (Int_t ilayer=0;ilayer<5;ilayer++){
1345       if (ntracks[ilayer]==0) continue;
1346       AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
1347       if (track.GetNumberOfClusters()<1) continue;
1348       CookLabel(&track,0);
1349       bestarray->AddAt(new AliITStrackMI(track),ilayer);
1350     }
1351   }
1352   //
1353   // update TPC V0 information
1354   //
1355   if (otrack->GetESDtrack()->GetV0Index(0)>0){
1356     Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
1357     for (Int_t i=0;i<3;i++){
1358       Int_t  index = otrack->GetESDtrack()->GetV0Index(i);
1359       if (index==0) break;
1360       AliV0 *vertex = (AliV0*)fEsd->GetV0(index);
1361       if (vertex->GetStatus()<0) continue;     // rejected V0
1362       //
1363       if (otrack->GetSign()>0) {
1364         vertex->SetIndex(0,esdindex);
1365       }
1366       else{
1367         vertex->SetIndex(1,esdindex);
1368       }
1369       //find nearest layer with track info
1370       Double_t xrp[3]; vertex->GetXYZ(xrp[0],xrp[1],xrp[2]);  //I.B.
1371       Int_t nearestold  = GetNearestLayer(xrp);               //I.B.
1372       Int_t nearest     = nearestold;
1373       for (Int_t ilayer =nearest;ilayer<8;ilayer++){
1374         if (ntracks[nearest]==0){
1375           nearest = ilayer;
1376         }
1377       }
1378       //
1379       AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
1380       if (nearestold<5&&nearest<5){
1381         Bool_t accept = track.GetNormChi2(nearest)<10;
1382         if (accept){
1383           if (track.GetSign()>0) {
1384             vertex->SetParamP(track);
1385             vertex->Update(fprimvertex);
1386             //vertex->SetIndex(0,track.fESDtrack->GetID());
1387             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1388           }else{
1389             vertex->SetParamN(track);
1390             vertex->Update(fprimvertex);
1391             //vertex->SetIndex(1,track.fESDtrack->GetID());
1392             if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
1393           }
1394           vertex->SetStatus(vertex->GetStatus()+1);
1395         }else{
1396           //vertex->SetStatus(-2);  // reject V0  - not enough clusters
1397         }
1398       }
1399     }
1400   }
1401
1402 }
1403 //------------------------------------------------------------------------
1404 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
1405 {
1406   //--------------------------------------------------------------------
1407   //
1408   //
1409   return fgLayers[layer];
1410 }
1411 //------------------------------------------------------------------------
1412 AliITStrackerMI::AliITSlayer::AliITSlayer():
1413 fR(0),
1414 fPhiOffset(0),
1415 fNladders(0),
1416 fZOffset(0),
1417 fNdetectors(0),
1418 fDetectors(0),
1419 fN(0),
1420 fDy5(0),
1421 fDy10(0),
1422 fDy20(0),
1423 fClustersCs(0),
1424 fClusterIndexCs(0),
1425 fYcs(0),
1426 fZcs(0),
1427 fNcs(0),
1428 fCurrentSlice(-1),
1429 fZmin(0),
1430 fZmax(0),
1431 fYmin(0),
1432 fYmax(0),
1433 fI(0),
1434 fImax(0),
1435 fSkip(0),
1436 fAccepted(0),
1437 fRoad(0),
1438 fMaxSigmaClY(0),
1439 fMaxSigmaClZ(0),
1440 fNMaxSigmaCl(3)
1441 {
1442   //--------------------------------------------------------------------
1443   //default AliITSlayer constructor
1444   //--------------------------------------------------------------------
1445   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1446     fClusterWeight[i]=0;
1447     fClusterTracks[0][i]=-1;
1448     fClusterTracks[1][i]=-1;
1449     fClusterTracks[2][i]=-1;
1450     fClusterTracks[3][i]=-1;
1451   }
1452 }
1453 //------------------------------------------------------------------------
1454 AliITStrackerMI::AliITSlayer::
1455 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd):
1456 fR(r),
1457 fPhiOffset(p),
1458 fNladders(nl),
1459 fZOffset(z),
1460 fNdetectors(nd),
1461 fDetectors(0),
1462 fN(0),
1463 fDy5(0),
1464 fDy10(0),
1465 fDy20(0),
1466 fClustersCs(0),
1467 fClusterIndexCs(0),
1468 fYcs(0),
1469 fZcs(0),
1470 fNcs(0),
1471 fCurrentSlice(-1),
1472 fZmin(0),
1473 fZmax(0),
1474 fYmin(0),
1475 fYmax(0),
1476 fI(0),
1477 fImax(0),
1478 fSkip(0),
1479 fAccepted(0),
1480 fRoad(0),
1481 fMaxSigmaClY(0),
1482 fMaxSigmaClZ(0),
1483 fNMaxSigmaCl(3) {
1484   //--------------------------------------------------------------------
1485   //main AliITSlayer constructor
1486   //--------------------------------------------------------------------
1487   fDetectors=new AliITSdetector[fNladders*fNdetectors];
1488   fRoad=2*fR*TMath::Sqrt(TMath::Pi()/1.);//assuming that there's only one cluster
1489 }
1490 //------------------------------------------------------------------------
1491 AliITStrackerMI::AliITSlayer::AliITSlayer(const AliITSlayer& layer):
1492 fR(layer.fR),
1493 fPhiOffset(layer.fPhiOffset),
1494 fNladders(layer.fNladders),
1495 fZOffset(layer.fZOffset),
1496 fNdetectors(layer.fNdetectors),
1497 fDetectors(layer.fDetectors),
1498 fN(layer.fN),
1499 fDy5(layer.fDy5),
1500 fDy10(layer.fDy10),
1501 fDy20(layer.fDy20),
1502 fClustersCs(layer.fClustersCs),
1503 fClusterIndexCs(layer.fClusterIndexCs),
1504 fYcs(layer.fYcs),
1505 fZcs(layer.fZcs),
1506 fNcs(layer.fNcs),
1507 fCurrentSlice(layer.fCurrentSlice),
1508 fZmin(layer.fZmin),
1509 fZmax(layer.fZmax),
1510 fYmin(layer.fYmin),
1511 fYmax(layer.fYmax),
1512 fI(layer.fI),
1513 fImax(layer.fImax),
1514 fSkip(layer.fSkip),
1515 fAccepted(layer.fAccepted),
1516 fRoad(layer.fRoad),
1517 fMaxSigmaClY(layer.fMaxSigmaClY),
1518 fMaxSigmaClZ(layer.fMaxSigmaClZ),
1519 fNMaxSigmaCl(layer.fNMaxSigmaCl)
1520 {
1521   //Copy constructor
1522 }
1523 //------------------------------------------------------------------------
1524 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
1525   //--------------------------------------------------------------------
1526   // AliITSlayer destructor
1527   //--------------------------------------------------------------------
1528   delete [] fDetectors;
1529   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1530   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1531     fClusterWeight[i]=0;
1532     fClusterTracks[0][i]=-1;
1533     fClusterTracks[1][i]=-1;
1534     fClusterTracks[2][i]=-1;
1535     fClusterTracks[3][i]=-1;
1536   }
1537 }
1538 //------------------------------------------------------------------------
1539 void AliITStrackerMI::AliITSlayer::ResetClusters() {
1540   //--------------------------------------------------------------------
1541   // This function removes loaded clusters
1542   //--------------------------------------------------------------------
1543   for (Int_t i=0; i<fN; i++) delete fClusters[i];
1544   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++){
1545     fClusterWeight[i]=0;
1546     fClusterTracks[0][i]=-1;
1547     fClusterTracks[1][i]=-1;
1548     fClusterTracks[2][i]=-1;
1549     fClusterTracks[3][i]=-1;
1550   }
1551
1552   fN=0;
1553   fI=0;
1554 }
1555 //------------------------------------------------------------------------
1556 void AliITStrackerMI::AliITSlayer::ResetWeights() {
1557   //--------------------------------------------------------------------
1558   // This function reset weights of the clusters
1559   //--------------------------------------------------------------------
1560   for (Int_t i=0; i<AliITSRecoParam::GetMaxClusterPerLayer(); i++) {
1561     fClusterWeight[i]=0;
1562     fClusterTracks[0][i]=-1;
1563     fClusterTracks[1][i]=-1;
1564     fClusterTracks[2][i]=-1;
1565     fClusterTracks[3][i]=-1;
1566   }
1567   for (Int_t i=0; i<fN;i++) {
1568     AliITSRecPoint * cl = (AliITSRecPoint*)GetCluster(i);
1569     if (cl&&cl->IsUsed()) cl->Use();
1570   }
1571
1572 }
1573 //------------------------------------------------------------------------
1574 void AliITStrackerMI::AliITSlayer::ResetRoad() {
1575   //--------------------------------------------------------------------
1576   // This function calculates the road defined by the cluster density
1577   //--------------------------------------------------------------------
1578   Int_t n=0;
1579   for (Int_t i=0; i<fN; i++) {
1580      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
1581   }
1582   if (n>1) fRoad=2*fR*TMath::Sqrt(TMath::Pi()/n);
1583 }
1584 //------------------------------------------------------------------------
1585 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSRecPoint *cl) {
1586   //--------------------------------------------------------------------
1587   //This function adds a cluster to this layer
1588   //--------------------------------------------------------------------
1589   if (fN==AliITSRecoParam::GetMaxClusterPerLayer()) {
1590     ::Error("InsertCluster","Too many clusters !\n");
1591     return 1;
1592   }
1593   fCurrentSlice=-1;
1594   fClusters[fN]=cl;
1595   fN++;
1596   AliITSdetector &det=GetDetector(cl->GetDetectorIndex());
1597   //AD
1598   Double_t nSigmaY=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaY2());
1599   Double_t nSigmaZ=fNMaxSigmaCl*TMath::Sqrt(cl->GetSigmaZ2());
1600   if (cl->GetY()-nSigmaY<det.GetYmin()) det.SetYmin(cl->GetY()-nSigmaY);
1601   if (cl->GetY()+nSigmaY>det.GetYmax()) det.SetYmax(cl->GetY()+nSigmaY);
1602   if (cl->GetZ()-nSigmaZ<det.GetZmin()) det.SetZmin(cl->GetZ()-nSigmaZ);
1603   if (cl->GetZ()+nSigmaZ>det.GetZmax()) det.SetZmax(cl->GetZ()+nSigmaZ);
1604   //AD
1605   /*
1606   if (cl->GetY()<det.GetYmin()) det.SetYmin(cl->GetY());
1607   if (cl->GetY()>det.GetYmax()) det.SetYmax(cl->GetY());
1608   if (cl->GetZ()<det.GetZmin()) det.SetZmin(cl->GetZ());
1609   if (cl->GetZ()>det.GetZmax()) det.SetZmax(cl->GetZ());
1610   */
1611   return 0;
1612 }
1613 //------------------------------------------------------------------------
1614 void  AliITStrackerMI::AliITSlayer::SortClusters()
1615 {
1616   //
1617   //sort clusters
1618   //
1619   AliITSRecPoint **clusters = new AliITSRecPoint*[fN];
1620   Float_t *z                = new Float_t[fN];
1621   Int_t   * index           = new Int_t[fN];
1622   //
1623   fMaxSigmaClY=0.; //AD
1624   fMaxSigmaClZ=0.; //AD
1625
1626   for (Int_t i=0;i<fN;i++){
1627     z[i] = fClusters[i]->GetZ();
1628     // save largest errors in y and z for this layer
1629     fMaxSigmaClY=TMath::Max(fMaxSigmaClY,TMath::Sqrt(fClusters[i]->GetSigmaY2()));
1630     fMaxSigmaClZ=TMath::Max(fMaxSigmaClZ,TMath::Sqrt(fClusters[i]->GetSigmaZ2()));
1631   }
1632   TMath::Sort(fN,z,index,kFALSE);
1633   for (Int_t i=0;i<fN;i++){
1634     clusters[i] = fClusters[index[i]];
1635   }
1636   //
1637   for (Int_t i=0;i<fN;i++){
1638     fClusters[i] = clusters[i];
1639     fZ[i]        = fClusters[i]->GetZ();
1640     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
1641     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
1642     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
1643     fY[i] = y;
1644   }
1645   delete[] index;
1646   delete[] z;
1647   delete[] clusters;
1648   //
1649
1650   fYB[0]=10000000;
1651   fYB[1]=-10000000;
1652   for (Int_t i=0;i<fN;i++){
1653     if (fY[i]<fYB[0]) fYB[0]=fY[i];
1654     if (fY[i]>fYB[1]) fYB[1]=fY[i];
1655     fClusterIndex[i] = i;
1656   }
1657   //
1658   // fill slices
1659   fDy5 = (fYB[1]-fYB[0])/5.;
1660   fDy10 = (fYB[1]-fYB[0])/10.;
1661   fDy20 = (fYB[1]-fYB[0])/20.;
1662   for (Int_t i=0;i<6;i++)  fN5[i] =0;
1663   for (Int_t i=0;i<11;i++) fN10[i]=0;
1664   for (Int_t i=0;i<21;i++) fN20[i]=0;
1665   //
1666   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;}
1667   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;}
1668   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;}
1669   //
1670   //
1671   for (Int_t i=0;i<fN;i++)
1672     for (Int_t irot=-1;irot<=1;irot++){
1673       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
1674       // slice 5
1675       for (Int_t slice=0; slice<6;slice++){
1676         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<AliITSRecoParam::GetMaxClusterPerLayer5()){
1677           fClusters5[slice][fN5[slice]] = fClusters[i];
1678           fY5[slice][fN5[slice]] = curY;
1679           fZ5[slice][fN5[slice]] = fZ[i];
1680           fClusterIndex5[slice][fN5[slice]]=i;
1681           fN5[slice]++;
1682         }
1683       }
1684       // slice 10
1685       for (Int_t slice=0; slice<11;slice++){
1686         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<AliITSRecoParam::GetMaxClusterPerLayer10()){
1687           fClusters10[slice][fN10[slice]] = fClusters[i];
1688           fY10[slice][fN10[slice]] = curY;
1689           fZ10[slice][fN10[slice]] = fZ[i];
1690           fClusterIndex10[slice][fN10[slice]]=i;
1691           fN10[slice]++;
1692         }
1693       }
1694       // slice 20
1695       for (Int_t slice=0; slice<21;slice++){
1696         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<AliITSRecoParam::GetMaxClusterPerLayer20()){
1697           fClusters20[slice][fN20[slice]] = fClusters[i];
1698           fY20[slice][fN20[slice]] = curY;
1699           fZ20[slice][fN20[slice]] = fZ[i];
1700           fClusterIndex20[slice][fN20[slice]]=i;
1701           fN20[slice]++;
1702         }
1703       }
1704     }
1705
1706   //
1707   // consistency check
1708   //
1709   for (Int_t i=0;i<fN-1;i++){
1710     if (fZ[i]>fZ[i+1]){
1711       printf("Bug\n");
1712     }
1713   }
1714   //
1715   for (Int_t slice=0;slice<21;slice++)
1716   for (Int_t i=0;i<fN20[slice]-1;i++){
1717     if (fZ20[slice][i]>fZ20[slice][i+1]){
1718       printf("Bug\n");
1719     }
1720   }
1721
1722
1723 }
1724 //------------------------------------------------------------------------
1725 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1726   //--------------------------------------------------------------------
1727   // This function returns the index of the nearest cluster
1728   //--------------------------------------------------------------------
1729   Int_t ncl=0;
1730   const Float_t *zcl;
1731   if (fCurrentSlice<0) {
1732     ncl = fN;
1733     zcl   = fZ;
1734   }
1735   else{
1736     ncl   = fNcs;
1737     zcl   = fZcs;;
1738   }
1739
1740   if (ncl==0) return 0;
1741   Int_t b=0, e=ncl-1, m=(b+e)/2;
1742   for (; b<e; m=(b+e)/2) {
1743     //    if (z > fClusters[m]->GetZ()) b=m+1;
1744     if (z > zcl[m]) b=m+1;
1745     else e=m;
1746   }
1747   return m;
1748 }
1749 //------------------------------------------------------------------------
1750 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 {
1751   //--------------------------------------------------------------------
1752   // This function computes the rectangular road for this track
1753   //--------------------------------------------------------------------
1754
1755
1756   AliITSdetector &det = fgLayers[ilayer].GetDetector(idet);
1757   // take into account the misalignment: propagate track to misaligned detector plane
1758   if (!track->Propagate(det.GetPhi(),det.GetRmisal())) return kFALSE;
1759
1760   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
1761                     TMath::Sqrt(track->GetSigmaZ2() +
1762                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1763                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
1764                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
1765   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
1766                     TMath::Sqrt(track->GetSigmaY2() +
1767                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1768                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
1769                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
1770
1771   // track at boundary between detectors, enlarge road
1772   Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidth();
1773   if ( (track->GetY()-dy < det.GetYmin()+boundaryWidth) ||
1774        (track->GetY()+dy > det.GetYmax()-boundaryWidth) ||
1775        (track->GetZ()-dz < det.GetZmin()+boundaryWidth) ||
1776        (track->GetZ()+dz > det.GetZmax()-boundaryWidth) ) {
1777     Float_t tgl = TMath::Abs(track->GetTgl());
1778     if (tgl > 1.) tgl=1.;
1779     Double_t deltaXNeighbDets=AliITSRecoParam::GetDeltaXNeighbDets();
1780     dz = TMath::Sqrt(dz*dz+deltaXNeighbDets*deltaXNeighbDets*tgl*tgl);
1781     Float_t snp = TMath::Abs(track->GetSnp());
1782     if (snp > AliITSReconstructor::GetRecoParam()->GetMaxSnp()) return kFALSE;
1783     dy = TMath::Sqrt(dy*dy+deltaXNeighbDets*deltaXNeighbDets*snp*snp);
1784   } // boundary
1785
1786   // add to the road a term (up to 2-3 mm) to deal with misalignments
1787   dy = TMath::Sqrt(dy*dy + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1788   dz = TMath::Sqrt(dz*dz + AliITSReconstructor::GetRecoParam()->GetRoadMisal()*AliITSReconstructor::GetRecoParam()->GetRoadMisal());
1789
1790   Double_t r = fgLayers[ilayer].GetR();
1791   zmin = track->GetZ() - dz;
1792   zmax = track->GetZ() + dz;
1793   ymin = track->GetY() + r*det.GetPhi() - dy;
1794   ymax = track->GetY() + r*det.GetPhi() + dy;
1795
1796   // bring track back to idead detector plane
1797   if (!track->Propagate(det.GetPhi(),det.GetR())) return kFALSE;
1798
1799   return kTRUE;
1800 }
1801 //------------------------------------------------------------------------
1802 void AliITStrackerMI::AliITSlayer::
1803 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1804   //--------------------------------------------------------------------
1805   // This function sets the "window"
1806   //--------------------------------------------------------------------
1807
1808   Double_t circle=2*TMath::Pi()*fR;
1809   fYmin = ymin;
1810   fYmax = ymax;
1811   fZmin = zmin;
1812   fZmax = zmax;
1813   // AD
1814   // enlarge road in y by maximum cluster error on this layer (3 sigma)
1815   fYmin -= fNMaxSigmaCl*fMaxSigmaClY;
1816   fYmax += fNMaxSigmaCl*fMaxSigmaClY;
1817   fZmin -= fNMaxSigmaCl*fMaxSigmaClZ;
1818   fZmax += fNMaxSigmaCl*fMaxSigmaClZ;
1819
1820   Float_t ymiddle = (fYmin+fYmax)*0.5;
1821   if (ymiddle<fYB[0]) {
1822     fYmin+=circle; fYmax+=circle; ymiddle+=circle;
1823   } else if (ymiddle>fYB[1]) {
1824     fYmin-=circle; fYmax-=circle; ymiddle-=circle;
1825   }
1826
1827   //
1828   fCurrentSlice =-1;
1829   // defualt take all
1830   fClustersCs = fClusters;
1831   fClusterIndexCs = fClusterIndex;
1832   fYcs  = fY;
1833   fZcs  = fZ;
1834   fNcs  = fN;
1835   //
1836   //is in 20 slice?
1837   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1838     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1839     if (slice<0) slice=0;
1840     if (slice>20) slice=20;
1841     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1842     if (isOK) {
1843       fCurrentSlice=slice;
1844       fClustersCs = fClusters20[fCurrentSlice];
1845       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1846       fYcs  = fY20[fCurrentSlice];
1847       fZcs  = fZ20[fCurrentSlice];
1848       fNcs  = fN20[fCurrentSlice];
1849     }
1850   }
1851   //
1852   //is in 10 slice?
1853   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1854     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1855     if (slice<0) slice=0;
1856     if (slice>10) slice=10;
1857     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1858     if (isOK) {
1859       fCurrentSlice=slice;
1860       fClustersCs = fClusters10[fCurrentSlice];
1861       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1862       fYcs  = fY10[fCurrentSlice];
1863       fZcs  = fZ10[fCurrentSlice];
1864       fNcs  = fN10[fCurrentSlice];
1865     }
1866   }
1867   //
1868   //is in 5 slice?
1869   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1870     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1871     if (slice<0) slice=0;
1872     if (slice>5) slice=5;
1873     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1874     if (isOK) {
1875       fCurrentSlice=slice;
1876       fClustersCs = fClusters5[fCurrentSlice];
1877       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1878       fYcs  = fY5[fCurrentSlice];
1879       fZcs  = fZ5[fCurrentSlice];
1880       fNcs  = fN5[fCurrentSlice];
1881     }
1882   }
1883   //
1884   fI        = FindClusterIndex(fZmin);
1885   fImax     = TMath::Min(FindClusterIndex(fZmax)+1,fNcs);
1886   fSkip     = 0;
1887   fAccepted = 0;
1888
1889   return;
1890 }
1891 //------------------------------------------------------------------------
1892 Int_t AliITStrackerMI::AliITSlayer::
1893 FindDetectorIndex(Double_t phi, Double_t z) const {
1894   //--------------------------------------------------------------------
1895   //This function finds the detector crossed by the track
1896   //--------------------------------------------------------------------
1897   Double_t dphi;
1898   if (fZOffset<0)            // old geometry
1899     dphi = -(phi-fPhiOffset);
1900   else                       // new geometry
1901     dphi = phi-fPhiOffset;
1902
1903
1904   if      (dphi <  0) dphi += 2*TMath::Pi();
1905   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1906   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1907   if (np>=fNladders) np-=fNladders;
1908   if (np<0)          np+=fNladders;
1909
1910
1911   Double_t dz=fZOffset-z;
1912   Double_t nnz = dz*(fNdetectors-1)*0.5/fZOffset+0.5;
1913   Int_t nz = (nnz<0 ? -1 : (Int_t)nnz);
1914   if (nz>=fNdetectors) return -1;
1915   if (nz<0)            return -1;
1916
1917   // ad hoc correction for 3rd ladder of SDD inner layer,
1918   // which is reversed (rotated by pi around local y)
1919   // this correction is OK only from AliITSv11Hybrid onwards
1920   if (GetR()>12. && GetR()<20.) { // SDD inner
1921     if(np==2) { // 3rd ladder
1922       nz = (fNdetectors-1) - nz;
1923     }
1924   }
1925   //printf("ndet %d phi %f z %f  np %d nz %d\n",fNdetectors,phi,z,np,nz);
1926
1927
1928   return np*fNdetectors + nz;
1929 }
1930 //------------------------------------------------------------------------
1931 const AliITSRecPoint *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci,Bool_t test)
1932 {
1933   //--------------------------------------------------------------------
1934   // This function returns clusters within the "window"
1935   //--------------------------------------------------------------------
1936
1937   if (fCurrentSlice<0) {
1938     Double_t rpi2 = 2.*fR*TMath::Pi();
1939     for (Int_t i=fI; i<fImax; i++) {
1940       Double_t y = fY[i];
1941       Double_t z = fZ[i];
1942       if (fYmax<y) y -= rpi2;
1943       if (fYmin>y) y += rpi2;
1944       if (y<fYmin) continue;
1945       if (y>fYmax) continue;
1946       // AD
1947       // skip clusters that are in "extended" road but they
1948       // 3sigma error does not touch the original road
1949       if (z+fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())<fZmin+fNMaxSigmaCl*fMaxSigmaClZ) continue;
1950       if (z-fNMaxSigmaCl*TMath::Sqrt(fClusters[i]->GetSigmaZ2())>fZmax-fNMaxSigmaCl*fMaxSigmaClZ) continue;
1951       //
1952       if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1953       ci=i;
1954       if (!test) fI=i+1;
1955       return fClusters[i];
1956     }
1957   } else {
1958     for (Int_t i=fI; i<fImax; i++) {
1959       if (fYcs[i]<fYmin) continue;
1960       if (fYcs[i]>fYmax) continue;
1961       if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1962       ci=fClusterIndexCs[i];
1963       if (!test) fI=i+1;
1964       return fClustersCs[i];
1965     }
1966   }
1967   return 0;
1968 }
1969 //------------------------------------------------------------------------
1970 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1971 const {
1972   //--------------------------------------------------------------------
1973   // This function returns the layer thickness at this point (units X0)
1974   //--------------------------------------------------------------------
1975   Double_t d=0.0085;
1976   x0=AliITSRecoParam::GetX0Air();
1977   if (43<fR&&fR<45) { //SSD2
1978      Double_t dd=0.0034;
1979      d=dd;
1980      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1981      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1982      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1983      for (Int_t i=0; i<12; i++) {
1984        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1985           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1986           d+=0.0034;
1987           break;
1988        }
1989        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1990           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1991           d+=0.0034;
1992           break;
1993        }
1994        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1995        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1996      }
1997   } else
1998   if (37<fR&&fR<41) { //SSD1
1999      Double_t dd=0.0034;
2000      d=dd;
2001      if (TMath::Abs(y-0.00)>3.40) d+=dd;
2002      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
2003      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
2004      for (Int_t i=0; i<11; i++) {
2005        if (TMath::Abs(z-3.9*i)<0.15) {
2006           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2007           d+=dd;
2008           break;
2009        }
2010        if (TMath::Abs(z+3.9*i)<0.15) {
2011           if (TMath::Abs(y-0.00)>3.40) d+=dd;
2012           d+=dd;
2013           break;
2014        }
2015        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2016        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
2017      }
2018   } else
2019   if (13<fR&&fR<26) { //SDD
2020      Double_t dd=0.0033;
2021      d=dd;
2022      if (TMath::Abs(y-0.00)>3.30) d+=dd;
2023
2024      if (TMath::Abs(y-1.80)<0.55) {
2025         d+=0.016;
2026         for (Int_t j=0; j<20; j++) {
2027           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2028           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2029         }
2030      }
2031      if (TMath::Abs(y+1.80)<0.55) {
2032         d+=0.016;
2033         for (Int_t j=0; j<20; j++) {
2034           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2035           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
2036         }
2037      }
2038
2039      for (Int_t i=0; i<4; i++) {
2040        if (TMath::Abs(z-7.3*i)<0.60) {
2041           d+=dd;
2042           if (TMath::Abs(y-0.00)>3.30) d+=dd;
2043           break;
2044        }
2045        if (TMath::Abs(z+7.3*i)<0.60) {
2046           d+=dd;
2047           if (TMath::Abs(y-0.00)>3.30) d+=dd;
2048           break;
2049        }
2050      }
2051   } else
2052   if (6<fR&&fR<8) {   //SPD2
2053      Double_t dd=0.0063; x0=21.5;
2054      d=dd;
2055      if (TMath::Abs(y-3.08)>0.5) d+=dd;
2056      if (TMath::Abs(y-3.03)<0.10) d+=0.014;
2057   } else
2058   if (3<fR&&fR<5) {   //SPD1
2059      Double_t dd=0.0063; x0=21.5;
2060      d=dd;
2061      if (TMath::Abs(y+0.21)>0.6) d+=dd;
2062      if (TMath::Abs(y+0.10)<0.10) d+=0.014;
2063   }
2064
2065   return d;
2066 }
2067 //------------------------------------------------------------------------
2068 AliITStrackerMI::AliITSdetector::AliITSdetector(const AliITSdetector& det):
2069 fR(det.fR),
2070 fRmisal(det.fRmisal),
2071 fPhi(det.fPhi),
2072 fSinPhi(det.fSinPhi),
2073 fCosPhi(det.fCosPhi),
2074 fYmin(det.fYmin),
2075 fYmax(det.fYmax),
2076 fZmin(det.fZmin),
2077 fZmax(det.fZmax),
2078 fIsBad(det.fIsBad),
2079 fNChips(det.fNChips),
2080 fChipIsBad(det.fChipIsBad)
2081 {
2082   //Copy constructor
2083 }
2084 //------------------------------------------------------------------------
2085 void AliITStrackerMI::AliITSdetector::ReadBadDetectorAndChips(Int_t ilayer,Int_t idet,
2086                                                const AliITSDetTypeRec *detTypeRec)
2087 {
2088   //--------------------------------------------------------------------
2089   // Read bad detectors and chips from calibration objects in AliITSDetTypeRec
2090   //--------------------------------------------------------------------
2091
2092   // In AliITSDetTypeRec, detector numbers go from 0 to 2197
2093   // while in the tracker they start from 0 for each layer
2094   for(Int_t il=0; il<ilayer; il++)
2095     idet += AliITSgeomTGeo::GetNLadders(il+1)*AliITSgeomTGeo::GetNDetectors(il+1);
2096
2097   Int_t detType;
2098   if (ilayer==0 || ilayer==1) {        // ----------  SPD
2099     detType = 0;
2100   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
2101     detType = 1;
2102   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
2103     detType = 2;
2104   } else {
2105     printf("AliITStrackerMI::AliITSdetector::InitBadFromOCDB: Wrong layer number %d\n",ilayer);
2106     return;
2107   }
2108
2109   // Get calibration from AliITSDetTypeRec
2110   AliITSCalibration *calib = (AliITSCalibration*)detTypeRec->GetCalibrationModel(idet);
2111   calib->SetModuleIndex(idet);
2112   AliITSCalibration *calibSPDdead = 0;
2113   if(detType==0) calibSPDdead = (AliITSCalibration*)detTypeRec->GetSPDDeadModel(idet); // TEMPORARY
2114   if (calib->IsBad() ||
2115       (detType==0 && calibSPDdead->IsBad())) // TEMPORARY
2116     {
2117       SetBad();
2118       //      printf("lay %d bad %d\n",ilayer,idet);
2119     }
2120
2121   // Get segmentation from AliITSDetTypeRec
2122   AliITSsegmentation *segm = (AliITSsegmentation*)detTypeRec->GetSegmentationModel(detType);
2123
2124   // Read info about bad chips
2125   fNChips = segm->GetMaximumChipIndex()+1;
2126   //printf("ilayer %d  detType %d idet %d fNChips %d %d  GetNumberOfChips %d\n",ilayer,detType,idet,fNChips,segm->GetMaximumChipIndex(),segm->GetNumberOfChips());
2127   if(fChipIsBad) { delete [] fChipIsBad; fChipIsBad=NULL; }
2128   fChipIsBad = new Bool_t[fNChips];
2129   for (Int_t iCh=0;iCh<fNChips;iCh++) {
2130     fChipIsBad[iCh] = calib->IsChipBad(iCh);
2131     if (detType==0 && calibSPDdead->IsChipBad(iCh)) fChipIsBad[iCh] = kTRUE; // TEMPORARY
2132     //if(fChipIsBad[iCh]) {printf("lay %d det %d bad chip %d\n",ilayer,idet,iCh);}
2133   }
2134
2135   return;
2136 }
2137 //------------------------------------------------------------------------
2138 Double_t AliITStrackerMI::GetEffectiveThickness()
2139 {
2140   //--------------------------------------------------------------------
2141   // Returns the thickness between the current layer and the vertex (units X0)
2142   //--------------------------------------------------------------------
2143
2144   if(fUseTGeo!=0) {
2145     if(fxOverX0Layer[0]<0) BuildMaterialLUT("Layers");
2146     if(fxOverX0Shield[0]<0) BuildMaterialLUT("Shields");
2147     if(fxOverX0Pipe<0) BuildMaterialLUT("Pipe");
2148   }
2149
2150   // beam pipe
2151   Double_t dPipe = (fUseTGeo==0 ? AliITSRecoParam::GetdPipe() : fxOverX0Pipe);
2152   Double_t d=dPipe*AliITSRecoParam::GetrPipe()*AliITSRecoParam::GetrPipe();
2153
2154   // layers
2155   Double_t x0=0;
2156   Double_t xn=fgLayers[fI].GetR();
2157   for (Int_t i=0; i<fI; i++) {
2158     Double_t xi=fgLayers[i].GetR();
2159     Double_t dLayer = (fUseTGeo==0 ? fgLayers[i].GetThickness(0,0,x0) : fxOverX0Layer[i]);
2160     d+=dLayer*xi*xi;
2161   }
2162
2163   // shields
2164   if (fI>1) {
2165     Double_t dshieldSPD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(0) : fxOverX0Shield[0]);
2166     d+=dshieldSPD*AliITSRecoParam::GetrInsideShield(0)*AliITSRecoParam::GetrInsideShield(0);
2167   }
2168   if (fI>3) {
2169     Double_t dshieldSDD = (fUseTGeo==0 ? AliITSRecoParam::Getdshield(1) : fxOverX0Shield[1]);
2170     d+=dshieldSDD*AliITSRecoParam::GetrInsideShield(1)*AliITSRecoParam::GetrInsideShield(1);
2171   }
2172   return d/(xn*xn);
2173 }
2174 //------------------------------------------------------------------------
2175 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
2176   //-------------------------------------------------------------------
2177   // This function returns number of clusters within the "window"
2178   //--------------------------------------------------------------------
2179   Int_t ncl=0;
2180   for (Int_t i=fI; i<fN; i++) {
2181     const AliITSRecPoint *c=fClusters[i];
2182     if (c->GetZ() > fZmax) break;
2183     if (c->IsUsed()) continue;
2184     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
2185     Double_t y=fR*det.GetPhi() + c->GetY();
2186
2187     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
2188     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
2189
2190     if (y<fYmin) continue;
2191     if (y>fYmax) continue;
2192     ncl++;
2193   }
2194   return ncl;
2195 }
2196 //------------------------------------------------------------------------
2197 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2198                                 const AliITStrackMI *clusters,Bool_t extra, Bool_t planeeff)
2199 {
2200   //--------------------------------------------------------------------
2201   // This function refits the track "track" at the position "x" using
2202   // the clusters from "clusters"
2203   // If "extra"==kTRUE,
2204   //    the clusters from overlapped modules get attached to "track"
2205   // If "planeff"==kTRUE,
2206   //    special approach for plane efficiency evaluation is applyed
2207   //--------------------------------------------------------------------
2208
2209   Int_t index[AliITSgeomTGeo::kNLayers];
2210   Int_t k;
2211   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2212   Int_t nc=clusters->GetNumberOfClusters();
2213   for (k=0; k<nc; k++) {
2214     Int_t idx=clusters->GetClusterIndex(k);
2215     Int_t ilayer=(idx&0xf0000000)>>28;
2216     index[ilayer]=idx;
2217   }
2218
2219   return RefitAt(xx,track,index,extra,planeeff); // call the method below
2220 }
2221 //------------------------------------------------------------------------
2222 Bool_t AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *track,
2223                                 const Int_t *clusters,Bool_t extra, Bool_t planeeff)
2224 {
2225   //--------------------------------------------------------------------
2226   // This function refits the track "track" at the position "x" using
2227   // the clusters from array
2228   // If "extra"==kTRUE,
2229   //    the clusters from overlapped modules get attached to "track"
2230   // If "planeff"==kTRUE,
2231   //    special approach for plane efficiency evaluation is applyed
2232   //--------------------------------------------------------------------
2233   Int_t index[AliITSgeomTGeo::kNLayers];
2234   Int_t k;
2235   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
2236   //
2237   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
2238     index[k]=clusters[k];
2239   }
2240
2241   // special for cosmics: check which the innermost layer crossed
2242   // by the track
2243   Int_t innermostlayer=5;
2244   Double_t drphi = TMath::Abs(track->GetD(0.,0.));
2245   for(innermostlayer=0; innermostlayer<AliITSgeomTGeo::GetNLayers(); innermostlayer++) {
2246     if(drphi < fgLayers[innermostlayer].GetR()) break;
2247   }
2248   AliDebug(2,Form(" drphi  %f  innermost %d",drphi,innermostlayer));
2249
2250   Int_t modstatus=1; // found
2251   Float_t xloc,zloc;
2252   Int_t from, to, step;
2253   if (xx > track->GetX()) {
2254       from=innermostlayer; to=AliITSgeomTGeo::GetNLayers();
2255       step=+1;
2256   } else {
2257       from=AliITSgeomTGeo::GetNLayers()-1; to=innermostlayer-1;
2258       step=-1;
2259   }
2260   TString dir = (step>0 ? "outward" : "inward");
2261
2262   for (Int_t ilayer = from; ilayer != to; ilayer += step) {
2263      AliITSlayer &layer=fgLayers[ilayer];
2264      Double_t r=layer.GetR();
2265
2266      if (step<0 && xx>r) break;
2267
2268      // material between SSD and SDD, SDD and SPD
2269      Double_t hI=ilayer-0.5*step;
2270      if (TMath::Abs(hI-3.5)<0.01) // SDDouter
2271        if(!CorrectForShieldMaterial(track,"SDD",dir)) return kFALSE;
2272      if (TMath::Abs(hI-1.5)<0.01) // SPDouter
2273        if(!CorrectForShieldMaterial(track,"SPD",dir)) return kFALSE;
2274
2275
2276      Double_t oldGlobXYZ[3];
2277      if (!track->GetXYZ(oldGlobXYZ)) return kFALSE;
2278
2279      // continue if we are already beyond this layer
2280      Double_t oldGlobR = TMath::Sqrt(oldGlobXYZ[0]*oldGlobXYZ[0]+oldGlobXYZ[1]*oldGlobXYZ[1]);
2281      if(step>0 && oldGlobR > r) continue; // going outward
2282      if(step<0 && oldGlobR < r) continue; // going inward
2283
2284      Double_t phi,z;
2285      if (!track->GetPhiZat(r,phi,z)) return kFALSE;
2286
2287      Int_t idet=layer.FindDetectorIndex(phi,z);
2288
2289      // check if we allow a prolongation without point for large-eta tracks
2290      Int_t skip = CheckSkipLayer(track,ilayer,idet);
2291      if (skip==2) {
2292        modstatus = 4; // out in z
2293        if(LocalModuleCoord(ilayer,idet,track,xloc,zloc)) { // local module coords
2294          track->SetModuleIndexInfo(ilayer,idet,modstatus,xloc,zloc);
2295        }
2296        // cross layer
2297        // apply correction for material of the current layer
2298        // add time if going outward
2299        CorrectForLayerMaterial(track,ilayer,oldGlobXYZ,dir);
2300        continue;
2301      }
2302
2303      if (idet<0) return kFALSE;
2304
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 void AliITStrackerMI::SetForceSkippingOfLayer() {
4277   //-----------------------------------------------------------------
4278   // Check if we are forced to skip layers
4279   // either we set to skip them in RecoParam
4280   // or they were off during data-taking
4281   //-----------------------------------------------------------------
4282
4283   const AliEventInfo *eventInfo = GetEventInfo();
4284
4285   for(Int_t l=0; l<AliITSgeomTGeo::kNLayers; l++) {
4286     fForceSkippingOfLayer[l] = 0;
4287     // check reco param
4288     if(AliITSReconstructor::GetRecoParam()->GetLayersToSkip(l)) fForceSkippingOfLayer[l] = 1;
4289     // check run info
4290
4291     if(eventInfo) {
4292       AliDebug(2,Form("GetEventInfo->GetTriggerCluster: %s",eventInfo->GetTriggerCluster()));
4293       if(l==0 || l==1)  {
4294         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSPD")) fForceSkippingOfLayer[l] = 1;
4295       } else if(l==2 || l==3) {
4296         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSDD")) fForceSkippingOfLayer[l] = 1;
4297       } else {
4298         if(!strstr(eventInfo->GetTriggerCluster(),"ITSSSD")) fForceSkippingOfLayer[l] = 1;
4299       }
4300     }
4301   }
4302   return;
4303 }
4304 //------------------------------------------------------------------------
4305 Int_t AliITStrackerMI::CheckSkipLayer(const AliITStrackMI *track,
4306                                       Int_t ilayer,Int_t idet) const {
4307   //-----------------------------------------------------------------
4308   // This method is used to decide whether to allow a prolongation
4309   // without clusters, because we want to skip the layer.
4310   // In this case the return value is > 0:
4311   // return 1: the user requested to skip a layer
4312   // return 2: track outside z acceptance
4313   //-----------------------------------------------------------------
4314
4315   if (ForceSkippingOfLayer(ilayer)) return 1;
4316
4317   Int_t innerLayCanSkip=0; // was 2, changed on 05.11.2009
4318
4319   if (idet<0 &&  // out in z
4320       ilayer>innerLayCanSkip &&
4321       AliITSReconstructor::GetRecoParam()->GetExtendedEtaAcceptance()) {
4322     // check if track will cross SPD outer layer
4323     Double_t phiAtSPD2,zAtSPD2;
4324     if (track->GetPhiZat(fgLayers[1].GetR(),phiAtSPD2,zAtSPD2)) {
4325       if (TMath::Abs(zAtSPD2)<2.*AliITSRecoParam::GetSPDdetzlength()) return 2;
4326     }
4327     return 2; // always allow skipping, changed on 05.11.2009
4328   }
4329
4330   return 0;
4331 }
4332 //------------------------------------------------------------------------
4333 Int_t AliITStrackerMI::CheckDeadZone(AliITStrackMI *track,
4334                                      Int_t ilayer,Int_t idet,
4335                                      Double_t dz,Double_t dy,
4336                                      Bool_t noClusters) const {
4337   //-----------------------------------------------------------------
4338   // This method is used to decide whether to allow a prolongation
4339   // without clusters, because there is a dead zone in the road.
4340   // In this case the return value is > 0:
4341   // return 1: dead zone at z=0,+-7cm in SPD
4342   // return 2: all road is "bad" (dead or noisy) from the OCDB
4343   // return 3: at least a chip is "bad" (dead or noisy) from the OCDB
4344   // return 4: at least a single channel is "bad" (dead or noisy) from the OCDB
4345   //-----------------------------------------------------------------
4346
4347   // check dead zones at z=0,+-7cm in the SPD
4348   if (ilayer<2 && !AliITSReconstructor::GetRecoParam()->GetAddVirtualClustersInDeadZone()) {
4349     Double_t zmindead[3]={fSPDdetzcentre[0] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4350                           fSPDdetzcentre[1] + 0.5*AliITSRecoParam::GetSPDdetzlength(),
4351                           fSPDdetzcentre[2] + 0.5*AliITSRecoParam::GetSPDdetzlength()};
4352     Double_t zmaxdead[3]={fSPDdetzcentre[1] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4353                           fSPDdetzcentre[2] - 0.5*AliITSRecoParam::GetSPDdetzlength(),
4354                           fSPDdetzcentre[3] - 0.5*AliITSRecoParam::GetSPDdetzlength()};
4355     for (Int_t i=0; i<3; i++)
4356       if (track->GetZ()-dz<zmaxdead[i] && track->GetZ()+dz>zmindead[i]) {
4357         AliDebug(2,Form("crack SPD %d",ilayer));
4358         return 1;
4359       }
4360   }
4361
4362   // check bad zones from OCDB
4363   if (!AliITSReconstructor::GetRecoParam()->GetUseBadZonesFromOCDB()) return 0;
4364
4365   if (idet<0) return 0;
4366
4367   AliITSdetector &det=fgLayers[ilayer].GetDetector(idet);
4368
4369   Int_t detType=-1;
4370   Float_t detSizeFactorX=0.0001,detSizeFactorZ=0.0001;
4371   if (ilayer==0 || ilayer==1) {        // ----------  SPD
4372     detType = 0;
4373   } else if (ilayer==2 || ilayer==3) { // ----------  SDD
4374     detType = 1;
4375     detSizeFactorX *= 2.;
4376   } else if (ilayer==4 || ilayer==5) { // ----------  SSD
4377     detType = 2;
4378   }
4379   AliITSsegmentation *segm = (AliITSsegmentation*)fkDetTypeRec->GetSegmentationModel(detType);
4380   if (detType==2) segm->SetLayer(ilayer+1);
4381   Float_t detSizeX = detSizeFactorX*segm->Dx();
4382   Float_t detSizeZ = detSizeFactorZ*segm->Dz();
4383
4384   // check if the road overlaps with bad chips
4385   Float_t xloc,zloc;
4386   LocalModuleCoord(ilayer,idet,track,xloc,zloc);
4387   Float_t zlocmin = zloc-dz;
4388   Float_t zlocmax = zloc+dz;
4389   Float_t xlocmin = xloc-dy;
4390   Float_t xlocmax = xloc+dy;
4391   Int_t chipsInRoad[100];
4392
4393   // check if road goes out of detector
4394   Bool_t touchNeighbourDet=kFALSE;
4395   if (TMath::Abs(xlocmin)>0.5*detSizeX) {xlocmin=-0.5*detSizeX; touchNeighbourDet=kTRUE;}
4396   if (TMath::Abs(xlocmax)>0.5*detSizeX) {xlocmax=+0.5*detSizeX; touchNeighbourDet=kTRUE;}
4397   if (TMath::Abs(zlocmin)>0.5*detSizeZ) {zlocmin=-0.5*detSizeZ; touchNeighbourDet=kTRUE;}
4398   if (TMath::Abs(zlocmax)>0.5*detSizeZ) {zlocmax=+0.5*detSizeZ; touchNeighbourDet=kTRUE;}
4399   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));
4400
4401   // check if this detector is bad
4402   if (det.IsBad()) {
4403     AliDebug(2,Form("lay %d  bad detector %d",ilayer,idet));
4404     if(!touchNeighbourDet) {
4405       return 2; // all detectors in road are bad
4406     } else {
4407       return 3; // at least one is bad
4408     }
4409   }
4410
4411   Int_t nChipsInRoad = segm->GetChipsInLocalWindow(chipsInRoad,zlocmin,zlocmax,xlocmin,xlocmax);
4412   AliDebug(2,Form("lay %d nChipsInRoad %d",ilayer,nChipsInRoad));
4413   if (!nChipsInRoad) return 0;
4414
4415   Bool_t anyBad=kFALSE,anyGood=kFALSE;
4416   for (Int_t iCh=0; iCh<nChipsInRoad; iCh++) {
4417     if (chipsInRoad[iCh]<0 || chipsInRoad[iCh]>det.GetNChips()-1) continue;
4418     AliDebug(2,Form("  chip %d bad %d",chipsInRoad[iCh],(Int_t)det.IsChipBad(chipsInRoad[iCh])));
4419     if (det.IsChipBad(chipsInRoad[iCh])) {
4420       anyBad=kTRUE;
4421     } else {
4422       anyGood=kTRUE;
4423     }
4424   }
4425
4426   if (!anyGood) {
4427     if(!touchNeighbourDet) {
4428       AliDebug(2,"all bad in road");
4429       return 2;  // all chips in road are bad
4430     } else {
4431       return 3; // at least a bad chip in road
4432     }
4433   }
4434
4435   if (anyBad) {
4436     AliDebug(2,"at least a bad in road");
4437     return 3; // at least a bad chip in road
4438   }
4439
4440
4441   if (!AliITSReconstructor::GetRecoParam()->GetUseSingleBadChannelsFromOCDB()
4442       || !noClusters) return 0;
4443
4444   // There are no clusters in road: check if there is at least
4445   // a bad SPD pixel or SDD anode or SSD strips on both sides
4446
4447   Int_t idetInITS=idet;
4448   for(Int_t l=0;l<ilayer;l++) idetInITS+=AliITSgeomTGeo::GetNLadders(l+1)*AliITSgeomTGeo::GetNDetectors(l+1);
4449
4450   if (fITSChannelStatus->AnyBadInRoad(idetInITS,zlocmin,zlocmax,xlocmin,xlocmax)) {
4451     AliDebug(2,Form("Bad channel in det %d of layer %d\n",idet,ilayer));
4452     return 4;
4453   }
4454   //if (fITSChannelStatus->FractionOfBadInRoad(idet,zlocmin,zlocmax,xlocmin,xlocmax) > AliITSReconstructor::GetRecoParam()->GetMinFractionOfBadInRoad()) return 3;
4455
4456   return 0;
4457 }
4458 //------------------------------------------------------------------------
4459 Bool_t AliITStrackerMI::LocalModuleCoord(Int_t ilayer,Int_t idet,
4460                                        const AliITStrackMI *track,
4461                                        Float_t &xloc,Float_t &zloc) const {
4462   //-----------------------------------------------------------------
4463   // Gives position of track in local module ref. frame
4464   //-----------------------------------------------------------------
4465
4466   xloc=0.;
4467   zloc=0.;
4468
4469   if(idet<0) return kTRUE; // track out of z acceptance of layer
4470
4471   Int_t ndet=AliITSgeomTGeo::GetNDetectors(ilayer+1); // layers from 1 to 6
4472
4473   Int_t lad = Int_t(idet/ndet) + 1;
4474
4475   Int_t det = idet - (lad-1)*ndet + 1;
4476
4477   Double_t xyzGlob[3],xyzLoc[3];
4478
4479   AliITSdetector &detector = fgLayers[ilayer].GetDetector(idet);
4480   // take into account the misalignment: xyz at real detector plane
4481   if(!track->GetXYZAt(detector.GetRmisal(),GetBz(),xyzGlob)) return kFALSE;
4482
4483   if(!AliITSgeomTGeo::GlobalToLocal(ilayer+1,lad,det,xyzGlob,xyzLoc)) return kFALSE;
4484
4485   xloc = (Float_t)xyzLoc[0];
4486   zloc = (Float_t)xyzLoc[2];
4487
4488   return kTRUE;
4489 }
4490 //------------------------------------------------------------------------
4491 Bool_t AliITStrackerMI::IsOKForPlaneEff(const AliITStrackMI* track, const Int_t *clusters, Int_t ilayer) const {
4492 //
4493 // Method to be optimized further:
4494 // Aim: decide whether a track can be used for PlaneEff evaluation
4495 //      the decision is taken based on the track quality at the layer under study
4496 //      no information on the clusters on this layer has to be used
4497 //      The criterium is to reject tracks at boundaries between basic block (e.g. SPD chip)
4498 //      the cut is done on number of sigmas from the boundaries
4499 //
4500 //  Input: Actual track, layer [0,5] under study
4501 //  Output: none
4502 //  Return: kTRUE if this is a good track
4503 //
4504 // it will apply a pre-selection to obtain good quality tracks.
4505 // Here also  you will have the possibility to put a control on the
4506 // impact point of the track on the basic block, in order to exclude border regions
4507 // this will be done by calling a proper method of the AliITSPlaneEff class.
4508 //
4509 // input: AliITStrackMI* track, ilayer= layer number [0,5]
4510 // return: Bool_t   -> kTRUE if usable track, kFALSE if not usable.
4511 //
4512   Int_t index[AliITSgeomTGeo::kNLayers];
4513   Int_t k;
4514   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) index[k]=-1;
4515   //
4516   for (k=0; k<AliITSgeomTGeo::GetNLayers(); k++) {
4517     index[k]=clusters[k];
4518   }
4519
4520   if(!fPlaneEff)
4521     {AliWarning("IsOKForPlaneEff: null pointer to AliITSPlaneEff"); return kFALSE;}
4522   AliITSlayer &layer=fgLayers[ilayer];
4523   Double_t r=layer.GetR();
4524   AliITStrackMI tmp(*track);
4525
4526 // require a minimal number of cluster in other layers and eventually clusters in closest layers
4527   Int_t ncl=0;
4528   for(Int_t lay=AliITSgeomTGeo::kNLayers-1;lay>ilayer;lay--) {
4529     AliDebug(2,Form("trak=%d  lay=%d  ; index=%d ESD label= %d",tmp.GetLabel(),lay,
4530                     tmp.GetClIndex(lay),((AliESDtrack*)tmp.GetESDtrack())->GetLabel())) ;
4531     if (tmp.GetClIndex(lay)>=0) ncl++;
4532   }
4533   Bool_t nextout = kFALSE;
4534   if(ilayer==AliITSgeomTGeo::kNLayers-1) nextout=kTRUE; // you are already on the outermost layer
4535   else nextout = ((tmp.GetClIndex(ilayer+1)>=0)? kTRUE : kFALSE );
4536   Bool_t nextin = kFALSE;
4537   if(ilayer==0) nextin=kTRUE; // you are already on the innermost layer
4538   else nextin = ((index[ilayer-1]>=0)? kTRUE : kFALSE );
4539   if(ncl<AliITSgeomTGeo::kNLayers-(ilayer+1)-AliITSReconstructor::GetRecoParam()->GetMaxMissingClustersPlaneEff())
4540      return kFALSE;
4541   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInOuterLayerPlaneEff() && !nextout)  return kFALSE;
4542   if(AliITSReconstructor::GetRecoParam()->GetRequireClusterInInnerLayerPlaneEff() && !nextin)   return kFALSE;
4543   if(tmp.Pt() < AliITSReconstructor::GetRecoParam()->GetMinPtPlaneEff()) return kFALSE;
4544  //  if(AliITSReconstructor::GetRecoParam()->GetOnlyConstraintPlaneEff()  && !tmp.GetConstrain()) return kFALSE;
4545
4546 // detector number
4547   Double_t phi,z;
4548   if (!tmp.GetPhiZat(r,phi,z)) return kFALSE;
4549   Int_t idet=layer.FindDetectorIndex(phi,z);
4550   if(idet<0) { AliInfo(Form("cannot find detector"));
4551     return kFALSE;}
4552
4553   // here check if it has good Chi Square.
4554
4555   //propagate to the intersection with the detector plane
4556   const AliITSdetector &det=layer.GetDetector(idet);
4557   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return kFALSE;
4558
4559   Float_t locx; //
4560   Float_t locz; //
4561   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return kFALSE;
4562   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4563   if(key>fPlaneEff->Nblock()) return kFALSE;
4564   Float_t blockXmn,blockXmx,blockZmn,blockZmx;
4565   if (!fPlaneEff->GetBlockBoundaries(key,blockXmn,blockXmx,blockZmn,blockZmx)) return kFALSE;
4566   //***************
4567   // DEFINITION OF SEARCH ROAD FOR accepting a track
4568   //
4569   //For the time being they are hard-wired, later on from AliITSRecoParam
4570   // Double_t nsigx=AliITSRecoParam::GetNSigXFarFromBoundary();
4571   // Double_t nsigz=AliITSRecoParam::GetNSigZFarFromBoundary();
4572   Double_t nsigz=4;
4573   Double_t nsigx=4;
4574   Double_t dx=nsigx*TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4575   Double_t dz=nsigz*TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4576                                                 // done for RecPoints
4577
4578   // exclude tracks at boundary between detectors
4579   //Double_t boundaryWidth=AliITSRecoParam::GetBoundaryWidthPlaneEff();
4580   Double_t boundaryWidth=0; // for the time being hard-wired, later on from AliITSRecoParam
4581   AliDebug(2,Form("Tracking: track impact x=%f, y=%f, z=%f",tmp.GetX(), tmp.GetY(), tmp.GetZ()));
4582   AliDebug(2,Form("Local:    track impact x=%f, z=%f",locx,locz));
4583   AliDebug(2,Form("Search Road. Tracking: dy=%f , dz=%f",dx,dz));
4584
4585   if ( (locx-dx < blockXmn+boundaryWidth) ||
4586        (locx+dx > blockXmx-boundaryWidth) ||
4587        (locz-dz < blockZmn+boundaryWidth) ||
4588        (locz+dz > blockZmx-boundaryWidth) ) return kFALSE;
4589   return kTRUE;
4590 }
4591 //------------------------------------------------------------------------
4592 void AliITStrackerMI::UseTrackForPlaneEff(const AliITStrackMI* track, Int_t ilayer) {
4593 //
4594 // This Method has to be optimized! For the time-being it uses the same criteria
4595 // as those used in the search of extra clusters for overlapping modules.
4596 //
4597 // Method Purpose: estabilish whether a track has produced a recpoint or not
4598 //                 in the layer under study (For Plane efficiency)
4599 //
4600 // inputs: AliITStrackMI* track  (pointer to a usable track)
4601 // outputs: none
4602 // side effects: update (by 1 count) the Plane Efficiency statistics of the basic block
4603 //               traversed by this very track. In details:
4604 //               - if a cluster can be associated to the track then call
4605 //                  AliITSPlaneEff::UpDatePlaneEff(key,kTRUE);
4606 //               - if not, the AliITSPlaneEff::UpDatePlaneEff(key,kFALSE) is called
4607 //
4608   if(!fPlaneEff)
4609     {AliWarning("UseTrackForPlaneEff: null pointer to AliITSPlaneEff"); return;}
4610   AliITSlayer &layer=fgLayers[ilayer];
4611   Double_t r=layer.GetR();
4612   AliITStrackMI tmp(*track);
4613
4614 // detector number
4615   Double_t phi,z;
4616   if (!tmp.GetPhiZat(r,phi,z)) return;
4617   Int_t idet=layer.FindDetectorIndex(phi,z);
4618
4619   if(idet<0) { AliInfo(Form("cannot find detector"));
4620     return;}
4621
4622
4623 //propagate to the intersection with the detector plane
4624   const AliITSdetector &det=layer.GetDetector(idet);
4625   if (!tmp.Propagate(det.GetPhi(),det.GetR())) return;
4626
4627
4628 //***************
4629 // DEFINITION OF SEARCH ROAD FOR CLUSTERS SELECTION
4630 //
4631   Double_t dz=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadZ()*
4632                     TMath::Sqrt(tmp.GetSigmaZ2() +
4633                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4634                     AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4635                     AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer));
4636   Double_t dy=AliITSReconstructor::GetRecoParam()->GetNSigmaRoadY()*
4637                     TMath::Sqrt(tmp.GetSigmaY2() +
4638                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4639                     AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4640                     AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer));
4641
4642 // road in global (rphi,z) [i.e. in tracking ref. system]
4643   Double_t zmin = tmp.GetZ() - dz;
4644   Double_t zmax = tmp.GetZ() + dz;
4645   Double_t ymin = tmp.GetY() + r*det.GetPhi() - dy;
4646   Double_t ymax = tmp.GetY() + r*det.GetPhi() + dy;
4647
4648 // select clusters in road
4649   layer.SelectClusters(zmin,zmax,ymin,ymax);
4650
4651 // Define criteria for track-cluster association
4652   Double_t msz = tmp.GetSigmaZ2() +
4653   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4654   AliITSReconstructor::GetRecoParam()->GetNSigmaZLayerForRoadZ()*
4655   AliITSReconstructor::GetRecoParam()->GetSigmaZ2(ilayer);
4656   Double_t msy = tmp.GetSigmaY2() +
4657   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4658   AliITSReconstructor::GetRecoParam()->GetNSigmaYLayerForRoadY()*
4659   AliITSReconstructor::GetRecoParam()->GetSigmaY2(ilayer);
4660   if (tmp.GetConstrain()) {
4661     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZC();
4662     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYC();
4663   }  else {
4664     msz *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadZNonC();
4665     msy *= AliITSReconstructor::GetRecoParam()->GetNSigma2RoadYNonC();
4666   }
4667   msz = 1./msz; // 1/RoadZ^2
4668   msy = 1./msy; // 1/RoadY^2
4669 //
4670
4671   const AliITSRecPoint *cl=0; Int_t clidx=-1, ci=-1;
4672   Int_t idetc=-1;
4673   Double_t chi2trkcl=1000.*AliITSReconstructor::GetRecoParam()->GetMaxChi2();
4674   //Double_t  tolerance=0.2;
4675   /*while ((cl=layer.GetNextCluster(clidx))!=0) {
4676     idetc = cl->GetDetectorIndex();
4677     if(idet!=idetc) continue;
4678     //Int_t ilay = cl->GetLayer();
4679
4680     if (TMath::Abs(tmp.GetZ() - cl->GetZ()) > tolerance) continue;
4681     if (TMath::Abs(tmp.GetY() - cl->GetY()) > tolerance) continue;
4682
4683     Double_t chi2=tmp.GetPredictedChi2(cl);
4684     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; }
4685   }*/
4686   Float_t locx; //
4687   Float_t locz; //
4688   if(!LocalModuleCoord(ilayer,idet,&tmp,locx,locz)) return;
4689 //
4690   AliDebug(2,Form("ilayer= %d, idet=%d, x= %f, z=%f",ilayer,idet,locx,locz));
4691   UInt_t key=fPlaneEff->GetKeyFromDetLocCoord(ilayer,idet,locx,locz);
4692   if(key>fPlaneEff->Nblock()) return;
4693   Bool_t found=kFALSE;
4694   //if (ci>=0) {
4695   Double_t chi2;
4696   while ((cl=layer.GetNextCluster(clidx))!=0) {
4697     idetc = cl->GetDetectorIndex();
4698     if(idet!=idetc) continue;
4699     // here real control to see whether the cluster can be associated to the track.
4700     // cluster not associated to track
4701     if ( (tmp.GetZ()-cl->GetZ())*(tmp.GetZ()-cl->GetZ())*msz +
4702          (tmp.GetY()-cl->GetY())*(tmp.GetY()-cl->GetY())*msy   > 1. ) continue;
4703     // calculate track-clusters chi2
4704     chi2 = GetPredictedChi2MI(&tmp,cl,ilayer); // note that this method change track tmp
4705                                                // in particular, the error associated to the cluster
4706     //Double_t chi2 = tmp.GetPredictedChi(cl); // this method does not change track tmp
4707     // chi2 cut
4708     if (chi2 > AliITSReconstructor::GetRecoParam()->GetMaxChi2s(ilayer)) continue;
4709     found=kTRUE;
4710     if (chi2<chi2trkcl) { chi2trkcl=chi2; ci=clidx; } // this just to trace which cluster is selected
4711    // track->SetExtraCluster(ilayer,(ilayer<<28)+ci);
4712    // track->SetExtraModule(ilayer,idetExtra);
4713   }
4714   if(!fPlaneEff->UpDatePlaneEff(found,key))
4715        AliWarning(Form("UseTrackForPlaneEff: cannot UpDate PlaneEff for key=%d",key));
4716   if(fPlaneEff->GetCreateHistos()&&  AliITSReconstructor::GetRecoParam()->GetHistoPlaneEff()) {
4717     Float_t tr[4]={99999.,99999.,9999.,9999.};    // initialize to high values
4718     Float_t clu[4]={-99999.,-99999.,9999.,9999.}; // (in some cases GetCov fails)
4719     Int_t cltype[2]={-999,-999};
4720
4721     tr[0]=locx;
4722     tr[1]=locz;
4723     tr[2]=TMath::Sqrt(tmp.GetSigmaY2());  // those are precisions in the tracking reference system
4724     tr[3]=TMath::Sqrt(tmp.GetSigmaZ2());  // Use it also for the module reference system, as it is
4725
4726     if (found){
4727       clu[0]=layer.GetCluster(ci)->GetDetLocalX();
4728       clu[1]=layer.GetCluster(ci)->GetDetLocalZ();
4729       cltype[0]=layer.GetCluster(ci)->GetNy();
4730       cltype[1]=layer.GetCluster(ci)->GetNz();
4731
4732      // Without the following 6 lines you would retrieve the nominal error of a cluster (e.g. for the SPD:
4733      //  X->50/sqrt(12)=14 micron   Z->450/sqrt(12)= 120 micron)
4734      // Within AliTrackerMI/AliTrackMI the error on the cluster is associated to the AliITStrackMI (fSigmaY,Z)
4735      // It is computed properly by calling the method
4736      // AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSRecPoint *cluster,Int_t layer)
4737      // T
4738      //Double_t x=0.5*(tmp.GetX()+layer.GetCluster(ci)->GetX()); // Take into account the mis-alignment
4739       //if (tmp.PropagateTo(x,0.,0.)) {
4740         chi2=GetPredictedChi2MI(&tmp,layer.GetCluster(ci),ilayer);
4741         AliCluster c(*layer.GetCluster(ci));
4742         c.SetSigmaY2(tmp.GetSigmaY(ilayer)*tmp.GetSigmaY(ilayer));
4743         c.SetSigmaZ2(tmp.GetSigmaZ(ilayer)*tmp.GetSigmaZ(ilayer));
4744         //if (layer.GetCluster(ci)->GetGlobalCov(cov))  // by using this, instead, you got nominal cluster errors
4745         clu[2]=TMath::Sqrt(c.GetSigmaY2());
4746         clu[3]=TMath::Sqrt(c.GetSigmaZ2());
4747       //}
4748     }
4749     fPlaneEff->FillHistos(key,found,tr,clu,cltype);
4750   }
4751 return;
4752 }