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