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