ITS/AliITStrackerMI.cxx

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