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