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