ITS/AliITStrackerMI.cxx

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