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