ITS/AliITStrackerMI.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-------------------------------------------------------------------------
  17 //               Implementation of the ITS tracker class
  18 //    It reads AliITSclusterV2 clusters and creates AliITStrackMI tracks
  19 //                   and fills with them the ESD
  20 //          Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
  21 //     dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
  22 //
  23 //-------------------------------------------------------------------------
  24 #include "AliITSrecoV2.h"
  25 #include <TTree.h>
  26 #include "AliITSgeom.h"
  27 #include "AliESD.h"
  28 #include "AliITSclusterV2.h"
  29 #include "AliITStrackerMI.h"
  30 #include "TMatrixD.h"
  31 #include "TFile.h"
  32 #include "TTree.h"
  33 #include "AliHelix.h"
  34 #include "AliESDV0MI.h"
  35 #include "AliLog.h"
  36 #include "TTreeStream.h"
  37
  38 ClassImp(AliITStrackerMI)
  39
  40
  41
  42 AliITStrackerMI::AliITSlayer AliITStrackerMI::fgLayers[kMaxLayer]; // ITS layers
  43
  44 AliITStrackerMI::AliITStrackerMI(const AliITSgeom *geom) : AliTracker() {
  45   //--------------------------------------------------------------------
  46   //This is the AliITStrackerMI constructor
  47   //--------------------------------------------------------------------
  48   fCoeficients = 0;
  49   fAfterV0     = kFALSE;
  50   AliITSgeom *g=(AliITSgeom*)geom;
  51   Float_t x,y,z;
  52   Int_t i;
  53   for (i=1; i<kMaxLayer+1; i++) {
  54     Int_t nlad=g->GetNladders(i);
  55     Int_t ndet=g->GetNdetectors(i);
  56
  57     g->GetTrans(i,1,1,x,y,z);
  58     Double_t r=TMath::Sqrt(x*x + y*y);
  59     Double_t poff=TMath::ATan2(y,x);
  60     Double_t zoff=z;
  61
  62     g->GetTrans(i,1,2,x,y,z);
  63     r += TMath::Sqrt(x*x + y*y);
  64     g->GetTrans(i,2,1,x,y,z);
  65     r += TMath::Sqrt(x*x + y*y);
  66     g->GetTrans(i,2,2,x,y,z);
  67     r += TMath::Sqrt(x*x + y*y);
  68     r*=0.25;
  69
  70     new (fgLayers+i-1) AliITSlayer(r,poff,zoff,nlad,ndet);
  71
  72     for (Int_t j=1; j<nlad+1; j++) {
  73       for (Int_t k=1; k<ndet+1; k++) { //Fill this layer with detectors
  74         Float_t x,y,zshift; g->GetTrans(i,j,k,x,y,zshift);
  75         Double_t rot[9]; g->GetRotMatrix(i,j,k,rot);
  76
  77         Double_t phi=TMath::ATan2(rot[1],rot[0])+TMath::Pi();
  78         phi+=TMath::Pi()/2;
  79         if (i==1) phi+=TMath::Pi();
  80         Double_t cp=TMath::Cos(phi), sp=TMath::Sin(phi);
  81         Double_t r=x*cp+y*sp;
  82
  83         AliITSdetector &det=fgLayers[i-1].GetDetector((j-1)*ndet + k-1);
  84         new(&det) AliITSdetector(r,phi);
  85       }
  86     }
  87
  88   }
  89
  90   fI=kMaxLayer;
  91
  92   fPass=0;
  93   fConstraint[0]=1; fConstraint[1]=0;
  94
  95   Double_t xyz[]={kXV,kYV,kZV}, ers[]={kSigmaXV,kSigmaYV,kSigmaZV};
  96   SetVertex(xyz,ers);
  97
  98   for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=kLayersNotToSkip[i];
  99   fLastLayerToTrackTo=kLastLayerToTrackTo;
 100   for (Int_t i=0;i<100000;i++){
 101     fBestTrackIndex[i]=0;
 102   }
 103 }
 104
 105 AliITStrackerMI::~AliITStrackerMI()
 106 {
 107   //
 108   //destructor
 109   //
 110   if (fCoeficients) delete []fCoeficients;
 111 }
 112
 113 void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
 114   //--------------------------------------------------------------------
 115   //This function set masks of the layers which must be not skipped
 116   //--------------------------------------------------------------------
 117   for (Int_t i=0; i<kMaxLayer; i++) fLayersNotToSkip[i]=l[i];
 118 }
 119
 120 Int_t AliITStrackerMI::LoadClusters(TTree *cTree) {
 121   //--------------------------------------------------------------------
 122   //This function loads ITS clusters
 123   //--------------------------------------------------------------------
 124   TBranch *branch=cTree->GetBranch("Clusters");
 125   if (!branch) {
 126     Error("LoadClusters"," can't get the branch !\n");
 127     return 1;
 128   }
 129
 130   TClonesArray dummy("AliITSclusterV2",10000), *clusters=&dummy;
 131   branch->SetAddress(&clusters);
 132
 133   Int_t j=0;
 134   Int_t detector=0;
 135   for (Int_t i=0; i<kMaxLayer; i++) {
 136     Int_t ndet=fgLayers[i].GetNdetectors();
 137     Int_t jmax = j + fgLayers[i].GetNladders()*ndet;
 138     for (; j<jmax; j++) {
 139       if (!cTree->GetEvent(j)) continue;
 140       Int_t ncl=clusters->GetEntriesFast();
 141       SignDeltas(clusters,GetZ());
 142       while (ncl--) {
 143         AliITSclusterV2 *c=(AliITSclusterV2*)clusters->UncheckedAt(ncl);
 144         detector = c->GetDetectorIndex();
 145         fgLayers[i].InsertCluster(new AliITSclusterV2(*c));
 146       }
 147       clusters->Delete();
 148       //add dead zone virtual "cluster"
 149       if (i<2){
 150         for (Float_t ydead = 0; ydead < 1.31 ; ydead+=(i+1.)*0.018){
 151           Int_t lab[4] = {0,0,0,detector};
 152           Int_t info[3] = {0,0,0};
 153           Float_t hit[5]={0,0,0.004/12.,0.001/12.,0};
 154           if (i==0) hit[0] =ydead-0.4;
 155           if (i==1) hit[0]=ydead-3.75;
 156           hit[1] =-0.04;
 157           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 158             fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
 159           hit[1]=-7.05;
 160           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 161             fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
 162           hit[1]=-7.15;
 163           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 164             fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
 165           hit[1] =0.06;
 166           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 167             fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
 168           hit[1]=7.05;
 169           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 170             fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
 171           hit[1]=7.25;
 172           if (TMath::Abs(fgLayers[i].GetDetector(detector).GetZmax()-hit[1])<2.)
 173             fgLayers[i].InsertCluster(new AliITSclusterV2(lab, hit, info));
 174         }
 175       }
 176
 177     }
 178     //
 179     fgLayers[i].ResetRoad(); //road defined by the cluster density
 180     fgLayers[i].SortClusters();
 181   }
 182
 183   return 0;
 184 }
 185
 186 void AliITStrackerMI::UnloadClusters() {
 187   //--------------------------------------------------------------------
 188   //This function unloads ITS clusters
 189   //--------------------------------------------------------------------
 190   for (Int_t i=0; i<kMaxLayer; i++) fgLayers[i].ResetClusters();
 191 }
 192
 193 static Int_t CorrectForDeadZoneMaterial(AliITStrackMI *t) {
 194   //--------------------------------------------------------------------
 195   // Correction for the material between the TPC and the ITS
 196   // (should it belong to the TPC code ?)
 197   //--------------------------------------------------------------------
 198   Double_t riw=80., diw=0.0053, x0iw=30; // TPC inner wall ?
 199   Double_t rcd=61., dcd=0.0053, x0cd=30; // TPC "central drum" ?
 200   Double_t yr=12.8, dr=0.03; // rods ?
 201   Double_t zm=0.2, dm=0.40;  // membrane
 202   //Double_t rr=52., dr=0.19, x0r=24., yyr=7.77; //rails
 203   Double_t rs=50., ds=0.001; // something belonging to the ITS (screen ?)
 204
 205   if (t->GetX() > riw) {
 206      if (!t->PropagateTo(riw,diw,x0iw)) return 1;
 207      if (TMath::Abs(t->GetY())>yr) t->CorrectForMaterial(dr);
 208      if (TMath::Abs(t->GetZ())<zm) t->CorrectForMaterial(dm);
 209      if (!t->PropagateTo(rcd,dcd,x0cd)) return 1;
 210      //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
 211      //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,dr,x0r);
 212      if (!t->PropagateTo(rs,ds)) return 1;
 213   } else if (t->GetX() < rs) {
 214      if (!t->PropagateTo(rs,-ds)) return 1;
 215      //Double_t x,y,z; t->GetGlobalXYZat(rr,x,y,z);
 216      //if (TMath::Abs(y)<yyr) t->PropagateTo(rr,-dr,x0r);
 217      if (!t->PropagateTo(rcd,-dcd,x0cd)) return 1;
 218      if (!t->PropagateTo(riw+0.001,-diw,x0iw)) return 1;
 219   } else {
 220   ::Error("CorrectForDeadZoneMaterial","track is already in the dead zone !");
 221     return 1;
 222   }
 223
 224   return 0;
 225 }
 226
 227 Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
 228   //--------------------------------------------------------------------
 229   // This functions reconstructs ITS tracks
 230   // The clusters must be already loaded !
 231   //--------------------------------------------------------------------
 232   TObjArray itsTracks(15000);
 233   fOriginal.Clear();
 234   {/* Read ESD tracks */
 235     Int_t nentr=event->GetNumberOfTracks();
 236     Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
 237     while (nentr--) {
 238       AliESDtrack *esd=event->GetTrack(nentr);
 239
 240       if ((esd->GetStatus()&AliESDtrack::kTPCin)==0) continue;
 241       if (esd->GetStatus()&AliESDtrack::kTPCout) continue;
 242       if (esd->GetStatus()&AliESDtrack::kITSin) continue;
 243       if (esd->GetKinkIndex(0)>0) continue;   //kink daughter
 244       AliITStrackMI *t=0;
 245       try {
 246         t=new AliITStrackMI(*esd);
 247       } catch (const Char_t *msg) {
 248         //Warning("Clusters2Tracks",msg);
 249         delete t;
 250         continue;
 251       }
 252       t->fD[0] = t->GetD(GetX(),GetY());
 253       t->fD[1] = t->GetZat(GetX())-GetZ();
 254       Double_t vdist = TMath::Sqrt(t->fD[0]*t->fD[0]+t->fD[1]*t->fD[1]);
 255       if (t->GetMass()<0.13) t->SetMass(0.13957); // MI look to the esd - mass hypothesys  !!!!!!!!!!!
 256       // write expected q
 257       t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.);
 258
 259       if (TMath::Abs(t->fD[0])>10) {
 260         delete t;
 261         continue;
 262       }
 263
 264       if (TMath::Abs(vdist)>20) {
 265         delete t;
 266         continue;
 267       }
 268       if (TMath::Abs(1/t->Get1Pt())<0.120) {
 269         delete t;
 270         continue;
 271       }
 272
 273       if (CorrectForDeadZoneMaterial(t)!=0) {
 274         //Warning("Clusters2Tracks",
 275         //        "failed to correct for the material in the dead zone !\n");
 276          delete t;
 277          continue;
 278       }
 279       t->fReconstructed = kFALSE;
 280       itsTracks.AddLast(t);
 281       fOriginal.AddLast(t);
 282     }
 283   } /* End Read ESD tracks */
 284
 285   itsTracks.Sort();
 286   fOriginal.Sort();
 287   Int_t nentr=itsTracks.GetEntriesFast();
 288   fTrackHypothesys.Expand(nentr);
 289   MakeCoeficients(nentr);
 290   Int_t ntrk=0;
 291   for (fPass=0; fPass<2; fPass++) {
 292      Int_t &constraint=fConstraint[fPass]; if (constraint<0) continue;
 293      for (Int_t i=0; i<nentr; i++) {
 294 //       cerr<<fPass<<"    "<<i<<'\r';
 295        fCurrentEsdTrack = i;
 296        AliITStrackMI *t=(AliITStrackMI*)itsTracks.UncheckedAt(i);
 297        if (t==0) continue;              //this track has been already tracked
 298        if (t->fReconstructed&&(t->fNUsed<1.5)) continue;  //this track was  already  "succesfully" reconstructed
 299        if ( (TMath::Abs(t->GetD(GetX(),GetY()))  >3.) && fConstraint[fPass]) continue;
 300        if ( (TMath::Abs(t->GetZat(GetX())-GetZ())>3.) && fConstraint[fPass]) continue;
 301
 302        Int_t tpcLabel=t->GetLabel(); //save the TPC track label
 303        fI = 6;
 304        ResetTrackToFollow(*t);
 305        ResetBestTrack();
 306
 307        FollowProlongationTree(t,i);
 308
 309
 310        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
 311        //
 312        AliITStrackMI * besttrack = GetBestHypothesys(fCurrentEsdTrack,t,15);
 313        if (!besttrack) continue;
 314        besttrack->SetLabel(tpcLabel);
 315        //       besttrack->CookdEdx();
 316        CookdEdx(besttrack);
 317        besttrack->fFakeRatio=1.;
 318        CookLabel(besttrack,0.); //For comparison only
 319        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
 320
 321        /*
 322        if ( besttrack->GetNumberOfClusters()<6 && fConstraint[fPass]) {
 323          continue;
 324        }
 325        if (besttrack->fChi2MIP[0]+besttrack->fNUsed>3.5) continue;
 326        if ( (TMath::Abs(besttrack->fD[0]*besttrack->fD[0]+besttrack->fD[1]*besttrack->fD[1])>0.1) && fConstraint[fPass])  continue;
 327        //delete itsTracks.RemoveAt(i);
 328        */
 329        if (fConstraint[fPass]&&(!besttrack->IsGoldPrimary())) continue;  //to be tracked also without vertex constrain
 330
 331
 332        t->fReconstructed = kTRUE;
 333        ntrk++;
 334      }
 335      GetBestHypothesysMIP(itsTracks);
 336   }
 337
 338   //GetBestHypothesysMIP(itsTracks);
 339   FindV0(event);
 340   fAfterV0 = kTRUE;
 341   //GetBestHypothesysMIP(itsTracks);
 342   //
 343   itsTracks.Delete();
 344   //
 345   Int_t entries = fTrackHypothesys.GetEntriesFast();
 346   for (Int_t ientry=0;ientry<entries;ientry++){
 347     TObjArray * array =(TObjArray*)fTrackHypothesys.UncheckedAt(ientry);
 348     if (array) array->Delete();
 349     delete fTrackHypothesys.RemoveAt(ientry);
 350   }
 351
 352   fTrackHypothesys.Delete();
 353   fOriginal.Clear();
 354   delete []fCoeficients;
 355   fCoeficients=0;
 356   Info("Clusters2Tracks","Number of prolonged tracks: %d\n",ntrk);
 357
 358   return 0;
 359 }
 360
 361
 362 Int_t AliITStrackerMI::PropagateBack(AliESD *event) {
 363   //--------------------------------------------------------------------
 364   // This functions propagates reconstructed ITS tracks back
 365   // The clusters must be loaded !
 366   //--------------------------------------------------------------------
 367   Int_t nentr=event->GetNumberOfTracks();
 368   Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 369
 370   Int_t ntrk=0;
 371   for (Int_t i=0; i<nentr; i++) {
 372      AliESDtrack *esd=event->GetTrack(i);
 373
 374      if ((esd->GetStatus()&AliESDtrack::kITSin)==0) continue;
 375      if (esd->GetStatus()&AliESDtrack::kITSout) continue;
 376
 377      AliITStrackMI *t=0;
 378      try {
 379         t=new AliITStrackMI(*esd);
 380      } catch (const Char_t *msg) {
 381        //Warning("PropagateBack",msg);
 382         delete t;
 383         continue;
 384      }
 385      t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.);
 386
 387      ResetTrackToFollow(*t);
 388
 389      // propagete to vertex [SR, GSI 17.02.2003]
 390      // Start Time measurement [SR, GSI 17.02.2003], corrected by I.Belikov
 391      if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {
 392        if (fTrackToFollow.PropagateToVertex()) {
 393           fTrackToFollow.StartTimeIntegral();
 394        }
 395        fTrackToFollow.PropagateTo(3.,-0.0028,65.19);
 396      }
 397
 398      fTrackToFollow.ResetCovariance(); fTrackToFollow.ResetClusters();
 399      if (RefitAt(49.,&fTrackToFollow,t)) {
 400         if (CorrectForDeadZoneMaterial(&fTrackToFollow)!=0) {
 401           //Warning("PropagateBack",
 402           //        "failed to correct for the material in the dead zone !\n");
 403           delete t;
 404           continue;
 405         }
 406         fTrackToFollow.SetLabel(t->GetLabel());
 407         //fTrackToFollow.CookdEdx();
 408         CookLabel(&fTrackToFollow,0.); //For comparison only
 409         fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSout);
 410         //UseClusters(&fTrackToFollow);
 411         ntrk++;
 412      }
 413      delete t;
 414   }
 415
 416   Info("PropagateBack","Number of back propagated ITS tracks: %d\n",ntrk);
 417
 418   return 0;
 419 }
 420
 421 Int_t AliITStrackerMI::RefitInward(AliESD *event) {
 422   //--------------------------------------------------------------------
 423   // This functions refits ITS tracks using the
 424   // "inward propagated" TPC tracks
 425   // The clusters must be loaded !
 426   //--------------------------------------------------------------------
 427   Int_t nentr=event->GetNumberOfTracks();
 428   Info("RefitInward", "Number of ESD tracks: %d\n", nentr);
 429
 430   Int_t ntrk=0;
 431   for (Int_t i=0; i<nentr; i++) {
 432     AliESDtrack *esd=event->GetTrack(i);
 433
 434     if ((esd->GetStatus()&AliESDtrack::kITSout) == 0) continue;
 435     if (esd->GetStatus()&AliESDtrack::kITSrefit) continue;
 436     if (esd->GetStatus()&AliESDtrack::kTPCout)
 437       if ((esd->GetStatus()&AliESDtrack::kTPCrefit)==0) continue;
 438
 439     AliITStrackMI *t=0;
 440     try {
 441         t=new AliITStrackMI(*esd);
 442     } catch (const Char_t *msg) {
 443       //Warning("RefitInward",msg);
 444         delete t;
 445         continue;
 446     }
 447     t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.);
 448     if (CorrectForDeadZoneMaterial(t)!=0) {
 449       //Warning("RefitInward",
 450       //         "failed to correct for the material in the dead zone !\n");
 451        delete t;
 452        continue;
 453     }
 454
 455     ResetTrackToFollow(*t);
 456     fTrackToFollow.ResetClusters();
 457
 458     if ((esd->GetStatus()&AliESDtrack::kTPCin)==0)
 459       fTrackToFollow.ResetCovariance();
 460
 461     //Refitting...
 462     if (RefitAt(3.7, &fTrackToFollow, t)) {
 463        fTrackToFollow.SetLabel(t->GetLabel());
 464        //       fTrackToFollow.CookdEdx();
 465        CookdEdx(&fTrackToFollow);
 466
 467        CookLabel(&fTrackToFollow,0.0); //For comparison only
 468
 469        if (fTrackToFollow.PropagateTo(3.,0.0028,65.19)) {//The beam pipe
 470          Double_t a=fTrackToFollow.GetAlpha();
 471          Double_t cs=TMath::Cos(a),sn=TMath::Sin(a);
 472          Double_t xv= GetX()*cs + GetY()*sn;
 473          Double_t yv=-GetX()*sn + GetY()*cs;
 474
 475          Double_t c=fTrackToFollow.GetC(), snp=fTrackToFollow.GetSnp();
 476          Double_t x=fTrackToFollow.GetX(), y=fTrackToFollow.GetY();
 477          Double_t tgfv=-(c*(x-xv)-snp)/(c*(y-yv) + TMath::Sqrt(1.-snp*snp));
 478          Double_t fv=TMath::ATan(tgfv);
 479
 480          cs=TMath::Cos(fv); sn=TMath::Sin(fv);
 481          x = xv*cs + yv*sn;
 482          yv=-xv*sn + yv*cs; xv=x;
 483
 484          if (fTrackToFollow.Propagate(fv+a,xv)) {
 485             fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
 486             Float_t d=fTrackToFollow.GetD(GetX(),GetY());
 487             Float_t z=fTrackToFollow.GetZ()-GetZ();
 488             fTrackToFollow.GetESDtrack()->SetImpactParameters(d,z);
 489             //UseClusters(&fTrackToFollow);
 490             {
 491             AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ());
 492             c.SetSigmaY2(GetSigmaY()*GetSigmaY());
 493             c.SetSigmaZ2(GetSigmaZ()*GetSigmaZ());
 494             Double_t chi2=fTrackToFollow.GetPredictedChi2(&c);
 495             //Double_t chi2=GetPredictedChi2MI(&fTrackToFollow,&c,fI);
 496             if (chi2<kMaxChi2)
 497               if (fTrackToFollow.Update(&c,-chi2,0))
 498                 //if (UpdateMI(&fTrackToFollow,&c,-chi2,0))
 499                 fTrackToFollow.SetConstrainedESDtrack(chi2);
 500             }
 501             ntrk++;
 502          }
 503        }
 504     }
 505     delete t;
 506   }
 507
 508   Info("RefitInward","Number of refitted tracks: %d\n",ntrk);
 509
 510   return 0;
 511 }
 512
 513 AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 514   //--------------------------------------------------------------------
 515   //       Return pointer to a given cluster
 516   //--------------------------------------------------------------------
 517   Int_t l=(index & 0xf0000000) >> 28;
 518   Int_t c=(index & 0x0fffffff) >> 00;
 519   return fgLayers[l].GetCluster(c);
 520 }
 521
 522
 523 void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex)
 524 {
 525   //--------------------------------------------------------------------
 526   // Follow prolongation tree
 527   //--------------------------------------------------------------------
 528
 529   //setup tree of the prolongations
 530   //
 531   static AliITStrackMI tracks[7][100];
 532   AliITStrackMI *currenttrack;
 533   static AliITStrackMI currenttrack1;
 534   static AliITStrackMI currenttrack2;
 535   static AliITStrackMI backuptrack;
 536   Int_t ntracks[7];
 537   Int_t nindexes[7][100];
 538   Float_t normalizedchi2[100];
 539   for (Int_t ilayer=0;ilayer<6;ilayer++) ntracks[ilayer]=0;
 540   otrack->fNSkipped=0;
 541   new (&(tracks[6][0])) AliITStrackMI(*otrack);
 542   ntracks[6]=1;
 543   nindexes[6][0]=0;
 544   //
 545   //
 546   // follow prolongations
 547   for (Int_t ilayer=5;ilayer>=0;ilayer--){
 548     //
 549     AliITSlayer &layer=fgLayers[ilayer];
 550     Double_t r=layer.GetR();
 551     ntracks[ilayer]=0;
 552     //
 553     //
 554    Int_t nskipped=0;
 555     Float_t nused =0;
 556     for (Int_t itrack =0;itrack<ntracks[ilayer+1];itrack++){
 557       //set current track
 558       if (ntracks[ilayer]>=100) break;
 559       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNSkipped>0) nskipped++;
 560       if (tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNUsed>2.) nused++;
 561       if (ntracks[ilayer]>15+ilayer){
 562         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNSkipped>0 && nskipped>4+ilayer) continue;
 563         if (itrack>1&&tracks[ilayer+1][nindexes[ilayer+1][itrack]].fNUsed>2. && nused>3) continue;
 564       }
 565
 566       new(&currenttrack1)  AliITStrackMI(tracks[ilayer+1][nindexes[ilayer+1][itrack]]);
 567       if (ilayer==3 || ilayer==1) {
 568         Double_t rs=0.5*(fgLayers[ilayer+1].GetR() + r);
 569         Double_t d=0.0034, x0=38.6;
 570         if (ilayer==1) {rs=9.; d=0.0097; x0=42;}
 571         if (!currenttrack1.PropagateTo(rs,d,x0)) {
 572           continue;
 573         }
 574       }
 575       //
 576       //find intersection with layer
 577       Double_t x,y,z;
 578       if (!currenttrack1.GetGlobalXYZat(r,x,y,z)) {
 579         continue;
 580       }
 581       Double_t phi=TMath::ATan2(y,x);
 582       Int_t idet=layer.FindDetectorIndex(phi,z);
 583       if (idet<0) {
 584         continue;
 585       }
 586       //propagate to the intersection
 587       const AliITSdetector &det=layer.GetDetector(idet);
 588       phi=det.GetPhi();
 589       new(&currenttrack2)  AliITStrackMI(currenttrack1);
 590       if (!currenttrack1.Propagate(phi,det.GetR())) {
 591         continue;
 592       }
 593       currenttrack2.Propagate(phi,det.GetR());  //
 594       currenttrack1.SetDetectorIndex(idet);
 595       currenttrack2.SetDetectorIndex(idet);
 596
 597       //
 598       //
 599       Double_t dz=7.5*TMath::Sqrt(currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]);
 600       Double_t dy=7.5*TMath::Sqrt(currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]);
 601       //
 602       Bool_t isBoundary=kFALSE;
 603       if (currenttrack1.GetY()-dy< det.GetYmin()+0.2) isBoundary = kTRUE;
 604       if (currenttrack1.GetY()+dy> det.GetYmax()-0.2) isBoundary = kTRUE;
 605       if (currenttrack1.GetZ()-dz< det.GetZmin()+0.2) isBoundary = kTRUE;
 606       if (currenttrack1.GetZ()+dz> det.GetZmax()-0.2) isBoundary = kTRUE;
 607
 608       if (isBoundary){ // track at boundary between detectors
 609         Float_t maxtgl = TMath::Abs(currenttrack1.GetTgl());
 610         if (maxtgl>1) maxtgl=1;
 611         dz = TMath::Sqrt(dz*dz+0.25*maxtgl*maxtgl);
 612         //
 613         Float_t maxsnp = TMath::Abs(currenttrack1.GetSnp());
 614         if (maxsnp>0.95) continue;
 615         //if (maxsnp>0.5) maxsnp=0.5;
 616         dy=TMath::Sqrt(dy*dy+0.25*maxsnp*maxsnp);
 617       }
 618
 619       Double_t zmin=currenttrack1.GetZ() - dz;
 620       Double_t zmax=currenttrack1.GetZ() + dz;
 621       Double_t ymin=currenttrack1.GetY() + r*phi - dy;
 622       Double_t ymax=currenttrack1.GetY() + r*phi + dy;
 623       layer.SelectClusters(zmin,zmax,ymin,ymax);
 624       //
 625       // loop over all possible prolongations
 626       //
 627       Double_t msz=1./((currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]));
 628       Double_t msy=1./((currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]));
 629       if (fConstraint[fPass]){
 630         msy/=60; msz/=60.;
 631       }
 632       else{
 633         msy/=50; msz/=50.;
 634       }
 635       //
 636       const AliITSclusterV2 *c=0; Int_t ci=-1;
 637       Double_t chi2=12345.;
 638       Int_t deadzone=0;
 639       currenttrack = &currenttrack1;
 640       while ((c=layer.GetNextCluster(ci))!=0) {
 641         if (ntracks[ilayer]>95) break; //space for skipped clusters
 642         Bool_t change =kFALSE;
 643         if (c->GetQ()==0 && (deadzone==1)) continue;
 644         Int_t idet=c->GetDetectorIndex();
 645         if (currenttrack->GetDetectorIndex()!=idet) {
 646           const AliITSdetector &det=layer.GetDetector(idet);
 647           Double_t y,z;
 648           if (!currenttrack2.GetProlongationFast(det.GetPhi(),det.GetR(),y,z)) continue;
 649           Float_t pz = (z - c->GetZ()) , py=(y - c->GetY());
 650           if (pz*pz*msz+py*py*msy>1.) continue;
 651           //
 652           new (&backuptrack) AliITStrackMI(currenttrack2);
 653           change = kTRUE;
 654           currenttrack =&currenttrack2;
 655           if (!currenttrack->Propagate(det.GetPhi(),det.GetR())) {
 656             new (currenttrack) AliITStrackMI(backuptrack);
 657             change = kFALSE;
 658             continue;
 659           }
 660           currenttrack->SetDetectorIndex(idet);
 661         }
 662         else{
 663           Float_t pz = (currenttrack->GetZ() - c->GetZ()) , py=(currenttrack->GetY() - c->GetY());
 664           if (pz*pz*msz+py*py*msy>1.) continue;
 665         }
 666
 667         chi2=GetPredictedChi2MI(currenttrack,c,ilayer);
 668         if (chi2<kMaxChi2s[ilayer]){
 669           if (c->GetQ()==0) deadzone=1;     // take dead zone only once
 670           if (ntracks[ilayer]>=100) continue;
 671           AliITStrackMI * updatetrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(*currenttrack);
 672           updatetrack->fClIndex[ilayer]=0;
 673           if (change){
 674             new (&currenttrack2) AliITStrackMI(backuptrack);
 675           }
 676           if (c->GetQ()!=0){
 677             if (!UpdateMI(updatetrack,c,chi2,(ilayer<<28)+ci)) continue;
 678             updatetrack->SetSampledEdx(c->GetQ(),updatetrack->GetNumberOfClusters()-1); //b.b.
 679           }
 680           else {
 681             updatetrack->fNDeadZone++;
 682             updatetrack->fDeadZoneProbability=GetDeadZoneProbability(updatetrack->GetZ(),TMath::Sqrt(updatetrack->GetSigmaZ2()));
 683           }
 684           if (c->IsUsed()){
 685             updatetrack->fNUsed++;
 686           }
 687           Double_t x0;
 688           Double_t d=layer.GetThickness(updatetrack->GetY(),updatetrack->GetZ(),x0);
 689           updatetrack->CorrectForMaterial(d,x0);
 690           if (fConstraint[fPass]) {
 691             updatetrack->fConstrain = fConstraint[fPass];
 692             fI = ilayer;
 693             Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 694             Double_t xyz[]={GetX(),GetY(),GetZ()};
 695             Double_t ptfactor = 1;
 696             Double_t ers[]={GetSigmaX()*ptfactor,GetSigmaY()*ptfactor,GetSigmaZ()};
 697             Bool_t isPrim = kTRUE;
 698             if (ilayer<4){
 699               updatetrack->fD[0] = updatetrack->GetD(GetX(),GetY());
 700               updatetrack->fD[1] = updatetrack->GetZat(GetX())-GetZ();
 701               if ( TMath::Abs(updatetrack->fD[0]/(1.+ilayer))>0.4 ||  TMath::Abs(updatetrack->fD[1]/(1.+ilayer))>0.4) isPrim=kFALSE;
 702             }
 703             if (isPrim) updatetrack->Improve(d,xyz,ers);
 704           } //apply vertex constrain
 705           ntracks[ilayer]++;
 706         }  // create new hypothesy
 707       } // loop over possible cluster prolongation
 708       //      if (fConstraint[fPass]&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0){
 709       if (fConstraint[fPass]&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){
 710         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
 711         vtrack->fClIndex[ilayer]=0;
 712         fI = ilayer;
 713         Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 714         Double_t xyz[]={GetX(),GetY(),GetZ()};
 715         Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 716         vtrack->Improve(d,xyz,ers);
 717         vtrack->fNSkipped++;
 718         ntracks[ilayer]++;
 719       }
 720
 721       if (fConstraint[fPass]&&itrack<1&&TMath::Abs(currenttrack1.fP3)>1.1){  //big theta -- for low mult. runs
 722         AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
 723         vtrack->fClIndex[ilayer]=0;
 724         fI = ilayer;
 725         Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
 726         Double_t xyz[]={GetX(),GetY(),GetZ()};
 727         Double_t ers[]={GetSigmaX(),GetSigmaY(),GetSigmaZ()};
 728         vtrack->Improve(d,xyz,ers);
 729         vtrack->fNDeadZone++;
 730         ntracks[ilayer]++;
 731       }
 732
 733
 734     } //loop over track candidates
 735     //
 736     //
 737     Int_t accepted=0;
 738
 739     Int_t golds=0;
 740     for (Int_t itrack=0;itrack<ntracks[ilayer];itrack++){
 741       normalizedchi2[itrack] = NormalizedChi2(&tracks[ilayer][itrack],ilayer);
 742       if ( normalizedchi2[itrack]<3+0.5*ilayer) golds++;
 743       if (ilayer>4) accepted++;
 744       else{
 745         if ( fConstraint[fPass] && normalizedchi2[itrack]<kMaxNormChi2C[ilayer]+1) accepted++;
 746         if (!fConstraint[fPass] && normalizedchi2[itrack]<kMaxNormChi2NonC[ilayer]+1) accepted++;
 747       }
 748     }
 749     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
 750     ntracks[ilayer] = TMath::Min(accepted,7+2*ilayer);
 751     if (ntracks[ilayer]<golds+2+ilayer) ntracks[ilayer]=TMath::Min(golds+2+ilayer,accepted);
 752     if (ntracks[ilayer]>90) ntracks[ilayer]=90;
 753   } //loop over layers
 754   //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]);
 755   Int_t max = fConstraint[fPass]? 20: 5;
 756
 757   for (Int_t i=0;i<TMath::Min(max,ntracks[0]);i++) {
 758     AliITStrackMI & track= tracks[0][nindexes[0][i]];
 759     if (track.GetNumberOfClusters()<2) continue;
 760     if (!fConstraint[fPass]&&track.fNormChi2[0]>7.)continue;
 761     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
 762   }
 763   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
 764     AliITStrackMI & track= tracks[1][nindexes[1][i]];
 765     if (track.GetNumberOfClusters()<4) continue;
 766     if (!fConstraint[fPass]&&track.fNormChi2[1]>7.)continue;
 767     if (fConstraint[fPass]) track.fNSkipped+=1;
 768     if (!fConstraint[fPass]) {
 769       track.fD[0] = track.GetD(GetX(),GetY());
 770       track.fNSkipped+=4./(4.+8.*TMath::Abs(track.fD[0]));
 771       if (track.fN+track.fNDeadZone+track.fNSkipped>6) {
 772         track.fNSkipped = 6-track.fN+track.fNDeadZone;
 773       }
 774     }
 775     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
 776   }
 777   //}
 778
 779   if (!fConstraint[fPass]){
 780     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
 781       AliITStrackMI & track= tracks[2][nindexes[2][i]];
 782       if (track.GetNumberOfClusters()<3) continue;
 783       if (!fConstraint[fPass]&&track.fNormChi2[2]>7.)continue;
 784       if (fConstraint[fPass]) track.fNSkipped+=2;
 785       if (!fConstraint[fPass]){
 786         track.fD[0] = track.GetD(GetX(),GetY());
 787         track.fNSkipped+= 7./(7.+8.*TMath::Abs(track.fD[0]));
 788         if (track.fN+track.fNDeadZone+track.fNSkipped>6) {
 789           track.fNSkipped = 6-track.fN+track.fNDeadZone;
 790         }
 791       }
 792       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
 793     }
 794   }
 795 }
 796
 797
 798 AliITStrackerMI::AliITSlayer & AliITStrackerMI::GetLayer(Int_t layer) const
 799 {
 800   //--------------------------------------------------------------------
 801   //
 802   //
 803   return fgLayers[layer];
 804 }
 805 AliITStrackerMI::AliITSlayer::AliITSlayer() {
 806   //--------------------------------------------------------------------
 807   //default AliITSlayer constructor
 808   //--------------------------------------------------------------------
 809   fN=0;
 810   fDetectors=0;
 811   fSkip = 0;
 812   fCurrentSlice=-1;
 813   for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
 814     fClusterWeight[i]=0;
 815     fClusterTracks[0][i]=-1;
 816     fClusterTracks[1][i]=-1;
 817     fClusterTracks[2][i]=-1;
 818     fClusterTracks[3][i]=-1;
 819   }
 820 }
 821
 822 AliITStrackerMI::AliITSlayer::
 823 AliITSlayer(Double_t r,Double_t p,Double_t z,Int_t nl,Int_t nd) {
 824   //--------------------------------------------------------------------
 825   //main AliITSlayer constructor
 826   //--------------------------------------------------------------------
 827   fR=r; fPhiOffset=p; fZOffset=z;
 828   fNladders=nl; fNdetectors=nd;
 829   fDetectors=new AliITSdetector[fNladders*fNdetectors];
 830
 831   fN=0;
 832   fI=0;
 833   fSkip = 0;
 834   fRoad=2*fR*TMath::Sqrt(3.14/1.);//assuming that there's only one cluster
 835 }
 836
 837 AliITStrackerMI::AliITSlayer::~AliITSlayer() {
 838   //--------------------------------------------------------------------
 839   // AliITSlayer destructor
 840   //--------------------------------------------------------------------
 841   delete[] fDetectors;
 842   for (Int_t i=0; i<fN; i++) delete fClusters[i];
 843   for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
 844     fClusterWeight[i]=0;
 845     fClusterTracks[0][i]=-1;
 846     fClusterTracks[1][i]=-1;
 847     fClusterTracks[2][i]=-1;
 848     fClusterTracks[3][i]=-1;
 849   }
 850 }
 851
 852 void AliITStrackerMI::AliITSlayer::ResetClusters() {
 853   //--------------------------------------------------------------------
 854   // This function removes loaded clusters
 855   //--------------------------------------------------------------------
 856   for (Int_t i=0; i<fN; i++) delete fClusters[i];
 857   for (Int_t i=0; i<kMaxClusterPerLayer;i++){
 858     fClusterWeight[i]=0;
 859     fClusterTracks[0][i]=-1;
 860     fClusterTracks[1][i]=-1;
 861     fClusterTracks[2][i]=-1;
 862     fClusterTracks[3][i]=-1;
 863   }
 864
 865   fN=0;
 866   fI=0;
 867 }
 868
 869 void AliITStrackerMI::AliITSlayer::ResetWeights() {
 870   //--------------------------------------------------------------------
 871   // This function reset weights of the clusters
 872   //--------------------------------------------------------------------
 873   for (Int_t i=0; i<kMaxClusterPerLayer;i++) {
 874     fClusterWeight[i]=0;
 875     fClusterTracks[0][i]=-1;
 876     fClusterTracks[1][i]=-1;
 877     fClusterTracks[2][i]=-1;
 878     fClusterTracks[3][i]=-1;
 879   }
 880   for (Int_t i=0; i<fN;i++) {
 881     AliITSclusterV2 * cl = (AliITSclusterV2*)GetCluster(i);
 882     if (cl&&cl->IsUsed()) cl->Use();
 883   }
 884
 885 }
 886
 887 void AliITStrackerMI::AliITSlayer::ResetRoad() {
 888   //--------------------------------------------------------------------
 889   // This function calculates the road defined by the cluster density
 890   //--------------------------------------------------------------------
 891   Int_t n=0;
 892   for (Int_t i=0; i<fN; i++) {
 893      if (TMath::Abs(fClusters[i]->GetZ())<fR) n++;
 894   }
 895   //if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
 896   if (n>1) fRoad=2*fR*TMath::Sqrt(3.14/n);
 897 }
 898
 899
 900 Int_t AliITStrackerMI::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
 901   //--------------------------------------------------------------------
 902   //This function adds a cluster to this layer
 903   //--------------------------------------------------------------------
 904   if (fN==kMaxClusterPerLayer) {
 905     ::Error("InsertCluster","Too many clusters !\n");
 906     return 1;
 907   }
 908   fCurrentSlice=-1;
 909   fClusters[fN]=c;
 910   fN++;
 911   AliITSdetector &det=GetDetector(c->GetDetectorIndex());
 912   if (c->GetY()<det.GetYmin()) det.SetYmin(c->GetY());
 913   if (c->GetY()>det.GetYmax()) det.SetYmax(c->GetY());
 914   if (c->GetZ()<det.GetZmin()) det.SetZmin(c->GetZ());
 915   if (c->GetZ()>det.GetZmax()) det.SetZmax(c->GetZ());
 916
 917   return 0;
 918 }
 919
 920 void  AliITStrackerMI::AliITSlayer::SortClusters()
 921 {
 922   //
 923   //sort clusters
 924   //
 925   AliITSclusterV2 **clusters = new AliITSclusterV2*[fN];
 926   Float_t *z                = new Float_t[fN];
 927   Int_t   * index           = new Int_t[fN];
 928   //
 929   for (Int_t i=0;i<fN;i++){
 930     z[i] = fClusters[i]->GetZ();
 931   }
 932   TMath::Sort(fN,z,index,kFALSE);
 933   for (Int_t i=0;i<fN;i++){
 934     clusters[i] = fClusters[index[i]];
 935   }
 936   //
 937   for (Int_t i=0;i<fN;i++){
 938     fClusters[i] = clusters[i];
 939     fZ[i]        = fClusters[i]->GetZ();
 940     AliITSdetector &det=GetDetector(fClusters[i]->GetDetectorIndex());
 941     Double_t y=fR*det.GetPhi() + fClusters[i]->GetY();
 942     if (y>2.*fR*TMath::Pi()) y -= 2.*fR*TMath::Pi();
 943     fY[i] = y;
 944   }
 945   delete[] index;
 946   delete[] z;
 947   delete[] clusters;
 948   //
 949
 950   fYB[0]=10000000;
 951   fYB[1]=-10000000;
 952   for (Int_t i=0;i<fN;i++){
 953     if (fY[i]<fYB[0]) fYB[0]=fY[i];
 954     if (fY[i]>fYB[1]) fYB[1]=fY[i];
 955     fClusterIndex[i] = i;
 956   }
 957   //
 958   // fill slices
 959   fDy5 = (fYB[1]-fYB[0])/5.;
 960   fDy10 = (fYB[1]-fYB[0])/10.;
 961   fDy20 = (fYB[1]-fYB[0])/20.;
 962   for (Int_t i=0;i<6;i++)  fN5[i] =0;
 963   for (Int_t i=0;i<11;i++) fN10[i]=0;
 964   for (Int_t i=0;i<21;i++) fN20[i]=0;
 965   //
 966   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;}
 967   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;}
 968   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;}
 969   //
 970   //
 971   for (Int_t i=0;i<fN;i++)
 972     for (Int_t irot=-1;irot<=1;irot++){
 973       Float_t curY = fY[i]+irot*TMath::TwoPi()*fR;
 974       // slice 5
 975       for (Int_t slice=0; slice<6;slice++){
 976         if (fBy5[slice][0]<curY && curY<fBy5[slice][1]&&fN5[slice]<kMaxClusterPerLayer5){
 977           fClusters5[slice][fN5[slice]] = fClusters[i];
 978           fY5[slice][fN5[slice]] = curY;
 979           fZ5[slice][fN5[slice]] = fZ[i];
 980           fClusterIndex5[slice][fN5[slice]]=i;
 981           fN5[slice]++;
 982         }
 983       }
 984       // slice 10
 985       for (Int_t slice=0; slice<11;slice++){
 986         if (fBy10[slice][0]<curY && curY<fBy10[slice][1]&&fN10[slice]<kMaxClusterPerLayer10){
 987           fClusters10[slice][fN10[slice]] = fClusters[i];
 988           fY10[slice][fN10[slice]] = curY;
 989           fZ10[slice][fN10[slice]] = fZ[i];
 990           fClusterIndex10[slice][fN10[slice]]=i;
 991           fN10[slice]++;
 992         }
 993       }
 994       // slice 20
 995       for (Int_t slice=0; slice<21;slice++){
 996         if (fBy20[slice][0]<curY && curY<fBy20[slice][1]&&fN20[slice]<kMaxClusterPerLayer20){
 997           fClusters20[slice][fN20[slice]] = fClusters[i];
 998           fY20[slice][fN20[slice]] = curY;
 999           fZ20[slice][fN20[slice]] = fZ[i];
1000           fClusterIndex20[slice][fN20[slice]]=i;
1001           fN20[slice]++;
1002         }
1003       }
1004     }
1005
1006   //
1007   // consistency check
1008   //
1009   for (Int_t i=0;i<fN-1;i++){
1010     if (fZ[i]>fZ[i+1]){
1011       printf("Bugg\n");
1012     }
1013   }
1014   //
1015   for (Int_t slice=0;slice<21;slice++)
1016   for (Int_t i=0;i<fN20[slice]-1;i++){
1017     if (fZ20[slice][i]>fZ20[slice][i+1]){
1018       printf("Bugg\n");
1019     }
1020   }
1021
1022
1023 }
1024
1025
1026 Int_t AliITStrackerMI::AliITSlayer::FindClusterIndex(Float_t z) const {
1027   //--------------------------------------------------------------------
1028   // This function returns the index of the nearest cluster
1029   //--------------------------------------------------------------------
1030   Int_t ncl=0;
1031   const Float_t *zcl;
1032   if (fCurrentSlice<0) {
1033     ncl = fN;
1034     zcl   = fZ;
1035   }
1036   else{
1037     ncl   = fNcs;
1038     zcl   = fZcs;;
1039   }
1040
1041   if (ncl==0) return 0;
1042   Int_t b=0, e=ncl-1, m=(b+e)/2;
1043   for (; b<e; m=(b+e)/2) {
1044     //    if (z > fClusters[m]->GetZ()) b=m+1;
1045     if (z > zcl[m]) b=m+1;
1046     else e=m;
1047   }
1048   return m;
1049 }
1050
1051
1052 void AliITStrackerMI::AliITSlayer::
1053 SelectClusters(Double_t zmin,Double_t zmax,Double_t ymin, Double_t ymax) {
1054   //--------------------------------------------------------------------
1055   // This function sets the "window"
1056   //--------------------------------------------------------------------
1057
1058   Double_t circle=2*TMath::Pi()*fR;
1059   fYmin = ymin; fYmax =ymax;
1060   Float_t ymiddle = (fYmin+fYmax)*0.5;
1061   if (ymiddle<fYB[0]) {fYmin+=circle; fYmax+=circle;ymiddle+=circle;}
1062   else{
1063     if (ymiddle>fYB[1]) {fYmin-=circle; fYmax-=circle;ymiddle-=circle;}
1064   }
1065   //
1066   fCurrentSlice =-1;
1067   // defualt take all
1068   fClustersCs = fClusters;
1069   fClusterIndexCs = fClusterIndex;
1070   fYcs  = fY;
1071   fZcs  = fZ;
1072   fNcs  = fN;
1073   //
1074   //is in 20 slice?
1075   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy20){
1076     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy20);
1077     if (slice<0) slice=0;
1078     if (slice>20) slice=20;
1079     Bool_t isOK = (fYmin>fBy20[slice][0]&&fYmax<fBy20[slice][1]);
1080     if (isOK) {
1081       fCurrentSlice=slice;
1082       fClustersCs = fClusters20[fCurrentSlice];
1083       fClusterIndexCs = fClusterIndex20[fCurrentSlice];
1084       fYcs  = fY20[fCurrentSlice];
1085       fZcs  = fZ20[fCurrentSlice];
1086       fNcs  = fN20[fCurrentSlice];
1087     }
1088   }
1089   //
1090   //is in 10 slice?
1091   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy10){
1092     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy10);
1093     if (slice<0) slice=0;
1094     if (slice>10) slice=10;
1095     Bool_t isOK = (fYmin>fBy10[slice][0]&&fYmax<fBy10[slice][1]);
1096     if (isOK) {
1097       fCurrentSlice=slice;
1098       fClustersCs = fClusters10[fCurrentSlice];
1099       fClusterIndexCs = fClusterIndex10[fCurrentSlice];
1100       fYcs  = fY10[fCurrentSlice];
1101       fZcs  = fZ10[fCurrentSlice];
1102       fNcs  = fN10[fCurrentSlice];
1103     }
1104   }
1105   //
1106   //is in 5 slice?
1107   if (fCurrentSlice<0&&TMath::Abs(fYmax-fYmin)<1.49*fDy5){
1108     Int_t slice = int(0.5+(ymiddle-fYB[0])/fDy5);
1109     if (slice<0) slice=0;
1110     if (slice>5) slice=5;
1111     Bool_t isOK = (fYmin>fBy5[slice][0]&&fYmax<fBy5[slice][1]);
1112     if ( isOK){
1113       fCurrentSlice=slice;
1114       fClustersCs = fClusters5[fCurrentSlice];
1115       fClusterIndexCs = fClusterIndex5[fCurrentSlice];
1116       fYcs  = fY5[fCurrentSlice];
1117       fZcs  = fZ5[fCurrentSlice];
1118       fNcs  = fN5[fCurrentSlice];
1119     }
1120   }
1121   //
1122   fI=FindClusterIndex(zmin); fZmax=zmax;
1123   fImax = TMath::Min(FindClusterIndex(zmax)+1,fNcs);
1124   fSkip = 0;
1125   fAccepted =0;
1126 }
1127
1128
1129
1130
1131 Int_t AliITStrackerMI::AliITSlayer::
1132 FindDetectorIndex(Double_t phi, Double_t z) const {
1133   //--------------------------------------------------------------------
1134   //This function finds the detector crossed by the track
1135   //--------------------------------------------------------------------
1136   Double_t dphi=-(phi-fPhiOffset);
1137   if      (dphi <  0) dphi += 2*TMath::Pi();
1138   else if (dphi >= 2*TMath::Pi()) dphi -= 2*TMath::Pi();
1139   Int_t np=Int_t(dphi*fNladders*0.5/TMath::Pi()+0.5);
1140   if (np>=fNladders) np-=fNladders;
1141   if (np<0)          np+=fNladders;
1142
1143   Double_t dz=fZOffset-z;
1144   Int_t nz=Int_t(dz*(fNdetectors-1)*0.5/fZOffset+0.5);
1145   if (nz>=fNdetectors) return -1;
1146   if (nz<0)            return -1;
1147
1148   return np*fNdetectors + nz;
1149 }
1150
1151
1152 const AliITSclusterV2 *AliITStrackerMI::AliITSlayer::GetNextCluster(Int_t &ci){
1153   //--------------------------------------------------------------------
1154   // This function returns clusters within the "window"
1155   //--------------------------------------------------------------------
1156
1157   if (fCurrentSlice<0){
1158     Double_t rpi2 = 2.*fR*TMath::Pi();
1159     for (Int_t i=fI; i<fImax; i++) {
1160       Double_t y = fY[i];
1161       if (fYmax<y) y -= rpi2;
1162       if (fYmin>y) y += rpi2;
1163       if (y<fYmin) continue;
1164       if (y>fYmax) continue;
1165       if (fClusters[i]->GetQ()==0&&fSkip==2) continue;
1166       ci=i;
1167       fI=i+1;
1168       return fClusters[i];
1169     }
1170   }
1171   else{
1172     for (Int_t i=fI; i<fImax; i++) {
1173       if (fYcs[i]<fYmin) continue;
1174       if (fYcs[i]>fYmax) continue;
1175       if (fClustersCs[i]->GetQ()==0&&fSkip==2) continue;
1176       ci=fClusterIndexCs[i];
1177       fI=i+1;
1178       return fClustersCs[i];
1179     }
1180   }
1181   return 0;
1182 }
1183
1184
1185
1186 Double_t AliITStrackerMI::AliITSlayer::GetThickness(Double_t y,Double_t z,Double_t &x0)
1187 const {
1188   //--------------------------------------------------------------------
1189   //This function returns the layer thickness at this point (units X0)
1190   //--------------------------------------------------------------------
1191   Double_t d=0.0085;
1192   x0=21.82;
1193   if (43<fR&&fR<45) { //SSD2
1194      Double_t dd=0.0034;
1195      d=dd;
1196      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1197      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1198      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1199      for (Int_t i=0; i<12; i++) {
1200        if (TMath::Abs(z-3.9*(i+0.5))<0.15) {
1201           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1202           d+=0.0034;
1203           break;
1204        }
1205        if (TMath::Abs(z+3.9*(i+0.5))<0.15) {
1206           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1207           d+=0.0034;
1208           break;
1209        }
1210        if (TMath::Abs(z-3.4-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1211        if (TMath::Abs(z+0.5+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1212      }
1213   } else
1214   if (37<fR&&fR<41) { //SSD1
1215      Double_t dd=0.0034;
1216      d=dd;
1217      if (TMath::Abs(y-0.00)>3.40) d+=dd;
1218      if (TMath::Abs(y-1.90)<0.45) {d+=(0.013-0.0034);}
1219      if (TMath::Abs(y+1.90)<0.45) {d+=(0.013-0.0034);}
1220      for (Int_t i=0; i<11; i++) {
1221        if (TMath::Abs(z-3.9*i)<0.15) {
1222           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1223           d+=dd;
1224           break;
1225        }
1226        if (TMath::Abs(z+3.9*i)<0.15) {
1227           if (TMath::Abs(y-0.00)>3.40) d+=dd;
1228           d+=dd;
1229           break;
1230        }
1231        if (TMath::Abs(z-1.85-3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1232        if (TMath::Abs(z+2.05+3.9*i)<0.50) {d+=(0.016-0.0034); break;}
1233      }
1234   } else
1235   if (13<fR&&fR<26) { //SDD
1236      Double_t dd=0.0033;
1237      d=dd;
1238      if (TMath::Abs(y-0.00)>3.30) d+=dd;
1239
1240      if (TMath::Abs(y-1.80)<0.55) {
1241         d+=0.016;
1242         for (Int_t j=0; j<20; j++) {
1243           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1244           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1245         }
1246      }
1247      if (TMath::Abs(y+1.80)<0.55) {
1248         d+=0.016;
1249         for (Int_t j=0; j<20; j++) {
1250           if (TMath::Abs(z-0.7-1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1251           if (TMath::Abs(z+0.7+1.47*j)<0.12) {d+=0.08; x0=9.; break;}
1252         }
1253      }
1254
1255      for (Int_t i=0; i<4; i++) {
1256        if (TMath::Abs(z-7.3*i)<0.60) {
1257           d+=dd;
1258           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1259           break;
1260        }
1261        if (TMath::Abs(z+7.3*i)<0.60) {
1262           d+=dd;
1263           if (TMath::Abs(y-0.00)>3.30) d+=dd;
1264           break;
1265        }
1266      }
1267   } else
1268   if (6<fR&&fR<8) {   //SPD2
1269      Double_t dd=0.0063; x0=21.5;
1270      d=dd;
1271      if (TMath::Abs(y-3.08)>0.5) d+=dd;
1272      //if (TMath::Abs(y-3.08)>0.45) d+=dd;
1273      if (TMath::Abs(y-3.03)<0.10) {d+=0.014;}
1274   } else
1275   if (3<fR&&fR<5) {   //SPD1
1276      Double_t dd=0.0063; x0=21.5;
1277      d=dd;
1278      if (TMath::Abs(y+0.21)>0.6) d+=dd;
1279      //if (TMath::Abs(y+0.21)>0.45) d+=dd;
1280      if (TMath::Abs(y+0.10)<0.10) {d+=0.014;}
1281   }
1282
1283   return d;
1284 }
1285
1286 Double_t AliITStrackerMI::GetEffectiveThickness(Double_t y,Double_t z) const
1287 {
1288   //--------------------------------------------------------------------
1289   //Returns the thickness between the current layer and the vertex (units X0)
1290   //--------------------------------------------------------------------
1291   Double_t d=0.0028*3*3; //beam pipe
1292   Double_t x0=0;
1293
1294   Double_t xn=fgLayers[fI].GetR();
1295   for (Int_t i=0; i<fI; i++) {
1296     Double_t xi=fgLayers[i].GetR();
1297     d+=fgLayers[i].GetThickness(y,z,x0)*xi*xi;
1298   }
1299
1300   if (fI>1) {
1301     Double_t xi=9.;
1302     d+=0.0097*xi*xi;
1303   }
1304
1305   if (fI>3) {
1306     Double_t xi=0.5*(fgLayers[3].GetR()+fgLayers[4].GetR());
1307     d+=0.0034*xi*xi;
1308   }
1309
1310   return d/(xn*xn);
1311 }
1312
1313 Int_t AliITStrackerMI::AliITSlayer::InRoad() const {
1314   //--------------------------------------------------------------------
1315   // This function returns number of clusters within the "window"
1316   //--------------------------------------------------------------------
1317   Int_t ncl=0;
1318   for (Int_t i=fI; i<fN; i++) {
1319     const AliITSclusterV2 *c=fClusters[i];
1320     if (c->GetZ() > fZmax) break;
1321     if (c->IsUsed()) continue;
1322     const AliITSdetector &det=GetDetector(c->GetDetectorIndex());
1323     Double_t y=fR*det.GetPhi() + c->GetY();
1324
1325     if (y>2.*fR*TMath::Pi()) y -= 2*fR*TMath::Pi();
1326     if (y>1.*fR*TMath::Pi() && fYmax<y) y -= 2*fR*TMath::Pi();
1327
1328     if (y<fYmin) continue;
1329     if (y>fYmax) continue;
1330     ncl++;
1331   }
1332   return ncl;
1333 }
1334
1335 Bool_t
1336 AliITStrackerMI::RefitAt(Double_t xx,AliITStrackMI *t,const AliITStrackMI *c) {
1337   //--------------------------------------------------------------------
1338   // This function refits the track "t" at the position "x" using
1339   // the clusters from "c"
1340   //--------------------------------------------------------------------
1341   Int_t index[kMaxLayer];
1342   Int_t k;
1343   for (k=0; k<kMaxLayer; k++) index[k]=-1;
1344   Int_t nc=c->GetNumberOfClusters();
1345   for (k=0; k<nc; k++) {
1346     Int_t idx=c->GetClusterIndex(k),nl=(idx&0xf0000000)>>28;
1347     index[nl]=idx;
1348   }
1349
1350   Int_t from, to, step;
1351   if (xx > t->GetX()) {
1352       from=0; to=kMaxLayer;
1353       step=+1;
1354   } else {
1355       from=kMaxLayer-1; to=-1;
1356       step=-1;
1357   }
1358
1359   for (Int_t i=from; i != to; i += step) {
1360      AliITSlayer &layer=fgLayers[i];
1361      Double_t r=layer.GetR();
1362
1363      {
1364      Double_t hI=i-0.5*step;
1365      if (TMath::Abs(hI-1.5)<0.01 || TMath::Abs(hI-3.5)<0.01) {
1366         Double_t rs=0.5*(fgLayers[i-step].GetR() + r);
1367         Double_t d=0.0034, x0=38.6;
1368         if (TMath::Abs(hI-1.5)<0.01) {rs=9.; d=0.0097; x0=42;}
1369         if (!t->PropagateTo(rs,-step*d,x0)) {
1370           return kFALSE;
1371         }
1372      }
1373      }
1374
1375      // remember old position [SR, GSI 18.02.2003]
1376      Double_t oldX=0., oldY=0., oldZ=0.;
1377      if (t->IsStartedTimeIntegral() && step==1) {
1378         t->GetGlobalXYZat(t->GetX(),oldX,oldY,oldZ);
1379      }
1380      //
1381
1382      Double_t x,y,z;
1383      if (!t->GetGlobalXYZat(r,x,y,z)) {
1384        return kFALSE;
1385      }
1386      Double_t phi=TMath::ATan2(y,x);
1387      Int_t idet=layer.FindDetectorIndex(phi,z);
1388      if (idet<0) {
1389        return kFALSE;
1390      }
1391      const AliITSdetector &det=layer.GetDetector(idet);
1392      phi=det.GetPhi();
1393      if (!t->Propagate(phi,det.GetR())) {
1394        return kFALSE;
1395      }
1396      t->SetDetectorIndex(idet);
1397
1398      const AliITSclusterV2 *cl=0;
1399      Double_t maxchi2=1000.*kMaxChi2;
1400
1401      Int_t idx=index[i];
1402      if (idx>0) {
1403         const AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(idx);
1404         if (idet != c->GetDetectorIndex()) {
1405            idet=c->GetDetectorIndex();
1406            const AliITSdetector &det=layer.GetDetector(idet);
1407            if (!t->Propagate(det.GetPhi(),det.GetR())) {
1408              return kFALSE;
1409            }
1410            t->SetDetectorIndex(idet);
1411         }
1412         //Double_t chi2=t->GetPredictedChi2(c);
1413         Int_t layer = (idx & 0xf0000000) >> 28;;
1414         Double_t chi2=GetPredictedChi2MI(t,c,layer);
1415         if (chi2<maxchi2) {
1416           cl=c;
1417           maxchi2=chi2;
1418         } else {
1419           return kFALSE;
1420         }
1421      }
1422      /*
1423      if (cl==0)
1424      if (t->GetNumberOfClusters()>2) {
1425         Double_t dz=4*TMath::Sqrt(t->GetSigmaZ2()+kSigmaZ2[i]);
1426         Double_t dy=4*TMath::Sqrt(t->GetSigmaY2()+kSigmaY2[i]);
1427         Double_t zmin=t->GetZ() - dz;
1428         Double_t zmax=t->GetZ() + dz;
1429         Double_t ymin=t->GetY() + phi*r - dy;
1430         Double_t ymax=t->GetY() + phi*r + dy;
1431         layer.SelectClusters(zmin,zmax,ymin,ymax);
1432
1433         const AliITSclusterV2 *c=0; Int_t ci=-1;
1434         while ((c=layer.GetNextCluster(ci))!=0) {
1435            if (idet != c->GetDetectorIndex()) continue;
1436            Double_t chi2=t->GetPredictedChi2(c);
1437            if (chi2<maxchi2) { cl=c; maxchi2=chi2; idx=ci; }
1438         }
1439      }
1440      */
1441      if (cl) {
1442        //if (!t->Update(cl,maxchi2,idx)) {
1443        if (!UpdateMI(t,cl,maxchi2,idx)) {
1444           return kFALSE;
1445        }
1446        t->SetSampledEdx(cl->GetQ(),t->GetNumberOfClusters()-1);
1447      }
1448
1449      {
1450      Double_t x0;
1451      Double_t d=layer.GetThickness(t->GetY(),t->GetZ(),x0);
1452      t->CorrectForMaterial(-step*d,x0);
1453      }
1454
1455      // track time update [SR, GSI 17.02.2003]
1456      if (t->IsStartedTimeIntegral() && step==1) {
1457         Double_t newX, newY, newZ;
1458         t->GetGlobalXYZat(t->GetX(),newX,newY,newZ);
1459         Double_t dL2 = (oldX-newX)*(oldX-newX) + (oldY-newY)*(oldY-newY) +
1460                        (oldZ-newZ)*(oldZ-newZ);
1461         t->AddTimeStep(TMath::Sqrt(dL2));
1462      }
1463      //
1464
1465   }
1466
1467   if (!t->PropagateTo(xx,0.,0.)) return kFALSE;
1468   return kTRUE;
1469 }
1470
1471
1472 Double_t AliITStrackerMI::GetNormalizedChi2(AliITStrackMI * track, Int_t mode)
1473 {
1474   //
1475   // calculate normalized chi2
1476   //  return NormalizedChi2(track,0);
1477   Float_t chi2 = 0;
1478   Float_t sum=0;
1479   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
1480   //  track->fdEdxMismatch=0;
1481   Float_t dedxmismatch =0;
1482   Float_t *ny = GetNy(fCurrentEsdTrack), *nz = GetNz(fCurrentEsdTrack);
1483   if (mode<100){
1484     for (Int_t i = 0;i<6;i++){
1485       if (track->fClIndex[i]>0){
1486         Float_t cerry, cerrz;
1487         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
1488         else
1489           { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
1490         cerry*=cerry;
1491         cerrz*=cerrz;
1492         Float_t cchi2 = (track->fDy[i]*track->fDy[i]/cerry)+(track->fDz[i]*track->fDz[i]/cerrz);
1493         if (i>1){
1494           Float_t ratio = track->fNormQ[i]/track->fExpQ;
1495           if (ratio<0.5) {
1496             cchi2+=(0.5-ratio)*10.;
1497             //track->fdEdxMismatch+=(0.5-ratio)*10.;
1498             dedxmismatch+=(0.5-ratio)*10.;
1499           }
1500         }
1501         if (i<2 ||i>3){
1502           AliITSclusterV2 * cl = (AliITSclusterV2*)GetCluster( track->fClIndex[i]);
1503           Double_t delta = cl->GetNy()+cl->GetNz()-ny[i]-nz[i];
1504           if (delta>1) chi2 +=0.5*TMath::Min(delta/2,2.);
1505           if (i<2) chi2+=2*cl->GetDeltaProbability();
1506         }
1507         chi2+=cchi2;
1508         sum++;
1509       }
1510     }
1511     if (TMath::Abs(track->fdEdxMismatch-dedxmismatch)>0.0001){
1512       track->fdEdxMismatch = dedxmismatch;
1513     }
1514   }
1515   else{
1516     for (Int_t i = 0;i<4;i++){
1517       if (track->fClIndex[i]>0){
1518         Float_t cerry, cerrz;
1519         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
1520         else { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
1521         cerry*=cerry;
1522         cerrz*=cerrz;
1523         chi2+= (track->fDy[i]*track->fDy[i]/cerry);
1524         chi2+= (track->fDz[i]*track->fDz[i]/cerrz);
1525         sum++;
1526       }
1527     }
1528     for (Int_t i = 4;i<6;i++){
1529       if (track->fClIndex[i]>0){
1530         Float_t cerry, cerrz;
1531         if (ny[i]>0) {cerry = erry[i]; cerrz=errz[i];}
1532         else { cerry= track->fSigmaY[i]; cerrz = track->fSigmaZ[i];}
1533         cerry*=cerry;
1534         cerrz*=cerrz;
1535         Float_t cerryb, cerrzb;
1536         if (ny[i+6]>0) {cerryb = erry[i+6]; cerrzb=errz[i+6];}
1537         else { cerryb= track->fSigmaY[i+6]; cerrzb = track->fSigmaZ[i+6];}
1538         cerryb*=cerryb;
1539         cerrzb*=cerrzb;
1540         chi2+= TMath::Min((track->fDy[i+6]*track->fDy[i+6]/cerryb),track->fDy[i]*track->fDy[i]/cerry);
1541         chi2+= TMath::Min((track->fDz[i+6]*track->fDz[i+6]/cerrzb),track->fDz[i]*track->fDz[i]/cerrz);
1542         sum++;
1543       }
1544     }
1545   }
1546   if (track->fESDtrack->GetTPCsignal()>85){
1547     Float_t ratio = track->fdEdx/track->fESDtrack->GetTPCsignal();
1548     if (ratio<0.5) {
1549       chi2+=(0.5-ratio)*5.;
1550     }
1551     if (ratio>2){
1552       chi2+=(ratio-2.0)*3;
1553     }
1554   }
1555   //
1556   Double_t match = TMath::Sqrt(track->fChi22);
1557   if (track->fConstrain)  match/=track->GetNumberOfClusters();
1558   if (!track->fConstrain) match/=track->GetNumberOfClusters()-2.;
1559   if (match<0) match=0;
1560   Float_t deadzonefactor = (track->fNDeadZone>0) ? 3*(1.1-track->fDeadZoneProbability):0.;
1561   Double_t normchi2 = 2*track->fNSkipped+match+deadzonefactor+(1+(2*track->fNSkipped+deadzonefactor)/track->GetNumberOfClusters())*
1562     (chi2)/TMath::Max(double(sum-track->fNSkipped),
1563                                 1./(1.+track->fNSkipped));
1564
1565  return normchi2;
1566 }
1567
1568
1569 Double_t AliITStrackerMI::GetMatchingChi2(AliITStrackMI * track1, AliITStrackMI * track2)
1570 {
1571   //
1572   // return matching chi2 between two tracks
1573   AliITStrackMI track3(*track2);
1574   track3.Propagate(track1->GetAlpha(),track1->GetX());
1575   TMatrixD vec(5,1);
1576   vec(0,0)=track1->fP0-track3.fP0;
1577   vec(1,0)=track1->fP1-track3.fP1;
1578   vec(2,0)=track1->fP2-track3.fP2;
1579   vec(3,0)=track1->fP3-track3.fP3;
1580   vec(4,0)=track1->fP4-track3.fP4;
1581   //
1582   TMatrixD cov(5,5);
1583   cov(0,0) = track1->fC00+track3.fC00;
1584   cov(1,1) = track1->fC11+track3.fC11;
1585   cov(2,2) = track1->fC22+track3.fC22;
1586   cov(3,3) = track1->fC33+track3.fC33;
1587   cov(4,4) = track1->fC44+track3.fC44;
1588
1589   cov(0,1)=cov(1,0) = track1->fC10+track3.fC10;
1590   cov(0,2)=cov(2,0) = track1->fC20+track3.fC20;
1591   cov(0,3)=cov(3,0) = track1->fC30+track3.fC30;
1592   cov(0,4)=cov(4,0) = track1->fC40+track3.fC40;
1593   //
1594   cov(1,2)=cov(2,1) = track1->fC21+track3.fC21;
1595   cov(1,3)=cov(3,1) = track1->fC31+track3.fC31;
1596   cov(1,4)=cov(4,1) = track1->fC41+track3.fC41;
1597   //
1598   cov(2,3)=cov(3,2) = track1->fC32+track3.fC32;
1599   cov(2,4)=cov(4,2) = track1->fC42+track3.fC42;
1600   //
1601   cov(3,4)=cov(4,3) = track1->fC43+track3.fC43;
1602
1603   cov.Invert();
1604   TMatrixD vec2(cov,TMatrixD::kMult,vec);
1605   TMatrixD chi2(vec2,TMatrixD::kTransposeMult,vec);
1606   return chi2(0,0);
1607 }
1608
1609 Double_t  AliITStrackerMI::GetDeadZoneProbability(Double_t zpos, Double_t zerr)
1610 {
1611   //
1612   //  return probability that given point - characterized by z position and error  is in dead zone
1613   //
1614   Double_t probability =0;
1615   Double_t absz = TMath::Abs(zpos);
1616   Double_t nearestz = (absz<2)? 0.:7.1;
1617   if (TMath::Abs(absz-nearestz)>0.25+3*zerr) return 0;
1618   Double_t zmin=0, zmax=0;
1619   if (zpos<-6.){
1620     zmin = -7.25; zmax = -6.95;
1621   }
1622   if (zpos>6){
1623     zmin = 7.0; zmax =7.3;
1624   }
1625   if (absz<2){
1626     zmin = -0.75; zmax = 1.5;
1627   }
1628   probability = (TMath::Erf((zpos-zmin)/zerr) - TMath::Erf((zpos-zmax)/zerr))*0.5;
1629   return probability;
1630 }
1631
1632
1633 Double_t AliITStrackerMI::GetTruncatedChi2(AliITStrackMI * track, Float_t fac)
1634 {
1635   //
1636   // calculate normalized chi2
1637   Float_t chi2[6];
1638   Float_t *erry = GetErrY(fCurrentEsdTrack), *errz = GetErrZ(fCurrentEsdTrack);
1639   Float_t ncl = 0;
1640   for (Int_t i = 0;i<6;i++){
1641     if (TMath::Abs(track->fDy[i])>0){
1642       chi2[i]= (track->fDy[i]/erry[i])*(track->fDy[i]/erry[i]);
1643       chi2[i]+= (track->fDz[i]/errz[i])*(track->fDz[i]/errz[i]);
1644       ncl++;
1645     }
1646     else{chi2[i]=10000;}
1647   }
1648   Int_t index[6];
1649   TMath::Sort(6,chi2,index,kFALSE);
1650   Float_t max = float(ncl)*fac-1.;
1651   Float_t sumchi=0, sumweight=0;
1652   for (Int_t i=0;i<max+1;i++){
1653     Float_t weight = (i<max)?1.:(max+1.-i);
1654     sumchi+=weight*chi2[index[i]];
1655     sumweight+=weight;
1656   }
1657   Double_t normchi2 = sumchi/sumweight;
1658   return normchi2;
1659 }
1660
1661
1662 Double_t AliITStrackerMI::GetInterpolatedChi2(AliITStrackMI * forwardtrack, AliITStrackMI * backtrack)
1663 {
1664   //
1665   // calculate normalized chi2
1666   //  if (forwardtrack->fNUsed>0.3*float(forwardtrack->GetNumberOfClusters())) return 10000;
1667   Int_t npoints = 0;
1668   Double_t res =0;
1669   for (Int_t i=0;i<6;i++){
1670     if ( (backtrack->fSigmaY[i]<0.000000001) || (forwardtrack->fSigmaY[i]<0.000000001)) continue;
1671     Double_t sy1 = forwardtrack->fSigmaY[i];
1672     Double_t sz1 = forwardtrack->fSigmaZ[i];
1673     Double_t sy2 = backtrack->fSigmaY[i];
1674     Double_t sz2 = backtrack->fSigmaZ[i];
1675     if (i<2){ sy2=1000.;sz2=1000;}
1676     //
1677     Double_t dy0 = (forwardtrack->fDy[i]/(sy1*sy1) +backtrack->fDy[i]/(sy2*sy2))/(1./(sy1*sy1)+1./(sy2*sy2));
1678     Double_t dz0 = (forwardtrack->fDz[i]/(sz1*sz1) +backtrack->fDz[i]/(sz2*sz2))/(1./(sz1*sz1)+1./(sz2*sz2));
1679     //
1680     Double_t nz0 = dz0*TMath::Sqrt((1./(sz1*sz1)+1./(sz2*sz2)));
1681     Double_t ny0 = dy0*TMath::Sqrt((1./(sy1*sy1)+1./(sy2*sy2)));
1682     //
1683     res+= nz0*nz0+ny0*ny0;
1684     npoints++;
1685   }
1686   if (npoints>1) return
1687                    TMath::Max(TMath::Abs(0.3*forwardtrack->Get1Pt())-0.5,0.)+
1688                    //2*forwardtrack->fNUsed+
1689                    res/TMath::Max(double(npoints-forwardtrack->fNSkipped),
1690                                   1./(1.+forwardtrack->fNSkipped));
1691   return 1000;
1692 }
1693
1694
1695
1696
1697
1698 Float_t  *AliITStrackerMI::GetWeight(Int_t index) {
1699   //--------------------------------------------------------------------
1700   //       Return pointer to a given cluster
1701   //--------------------------------------------------------------------
1702   Int_t l=(index & 0xf0000000) >> 28;
1703   Int_t c=(index & 0x0fffffff) >> 00;
1704   return fgLayers[l].GetWeight(c);
1705 }
1706
1707 void AliITStrackerMI::RegisterClusterTracks(AliITStrackMI* track,Int_t id)
1708 {
1709   //---------------------------------------------
1710   // register track to the list
1711   //
1712   if (track->fESDtrack->GetKinkIndex(0)!=0) return;  //don't register kink tracks
1713   //
1714   //
1715   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1716     Int_t index = track->GetClusterIndex(icluster);
1717     Int_t l=(index & 0xf0000000) >> 28;
1718     Int_t c=(index & 0x0fffffff) >> 00;
1719     if (c>fgLayers[l].fN) continue;
1720     for (Int_t itrack=0;itrack<4;itrack++){
1721       if (fgLayers[l].fClusterTracks[itrack][c]<0){
1722         fgLayers[l].fClusterTracks[itrack][c]=id;
1723         break;
1724       }
1725     }
1726   }
1727 }
1728 void AliITStrackerMI::UnRegisterClusterTracks(AliITStrackMI* track, Int_t id)
1729 {
1730   //---------------------------------------------
1731   // unregister track from the list
1732   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1733     Int_t index = track->GetClusterIndex(icluster);
1734     Int_t l=(index & 0xf0000000) >> 28;
1735     Int_t c=(index & 0x0fffffff) >> 00;
1736     if (c>fgLayers[l].fN) continue;
1737     for (Int_t itrack=0;itrack<4;itrack++){
1738       if (fgLayers[l].fClusterTracks[itrack][c]==id){
1739         fgLayers[l].fClusterTracks[itrack][c]=-1;
1740       }
1741     }
1742   }
1743 }
1744 Float_t AliITStrackerMI::GetNumberOfSharedClusters(AliITStrackMI* track,Int_t id, Int_t list[6], AliITSclusterV2 *clist[6])
1745 {
1746   //-------------------------------------------------------------
1747   //get number of shared clusters
1748   //-------------------------------------------------------------
1749   Float_t shared=0;
1750   for (Int_t i=0;i<6;i++) { list[i]=-1, clist[i]=0;}
1751   // mean  number of clusters
1752   Float_t *ny = GetNy(id), *nz = GetNz(id);
1753
1754
1755   for (Int_t icluster=0;icluster<track->GetNumberOfClusters();icluster++){
1756     Int_t index = track->GetClusterIndex(icluster);
1757     Int_t l=(index & 0xf0000000) >> 28;
1758     Int_t c=(index & 0x0fffffff) >> 00;
1759     if (c>fgLayers[l].fN) continue;
1760     if (ny[l]==0){
1761       printf("problem\n");
1762     }
1763     AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(index);
1764     Float_t weight=1;
1765     //
1766     Float_t deltan = 0;
1767     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
1768     if (l>2&&track->fNormQ[l]/track->fExpQ>3.5) continue;
1769     if (l<2 || l>3){
1770       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
1771     }
1772     else{
1773       deltan = (cl->GetNz()-nz[l]);
1774     }
1775     if (deltan>2.0) continue;  // extended - highly probable shared cluster
1776     weight = 2./TMath::Max(3.+deltan,2.);
1777     //
1778     for (Int_t itrack=0;itrack<4;itrack++){
1779       if (fgLayers[l].fClusterTracks[itrack][c]>=0 && fgLayers[l].fClusterTracks[itrack][c]!=id){
1780         list[l]=index;
1781         clist[l] = (AliITSclusterV2*)GetCluster(index);
1782         shared+=weight;
1783         break;
1784       }
1785     }
1786   }
1787   track->fNUsed=shared;
1788   return shared;
1789 }
1790
1791 Int_t AliITStrackerMI::GetOverlapTrack(AliITStrackMI *track, Int_t trackID, Int_t &shared, Int_t clusterlist[6],Int_t overlist[6])
1792 {
1793   //
1794   // find first shared track
1795   //
1796   // mean  number of clusters
1797   Float_t *ny = GetNy(trackID), *nz = GetNz(trackID);
1798   //
1799   for (Int_t i=0;i<6;i++) overlist[i]=-1;
1800   Int_t sharedtrack=100000;
1801   Int_t tracks[24],trackindex=0;
1802   for (Int_t i=0;i<24;i++) {tracks[i]=-1;}
1803   //
1804   for (Int_t icluster=0;icluster<6;icluster++){
1805     if (clusterlist[icluster]<0) continue;
1806     Int_t index = clusterlist[icluster];
1807     Int_t l=(index & 0xf0000000) >> 28;
1808     Int_t c=(index & 0x0fffffff) >> 00;
1809     if (ny[l]==0){
1810       printf("problem\n");
1811     }
1812     if (c>fgLayers[l].fN) continue;
1813     //if (l>3) continue;
1814     AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(index);
1815     //
1816     Float_t deltan = 0;
1817     if (l>3&&cl->GetNy()+cl->GetNz()>6) continue;
1818     if (l>2&&track->fNormQ[l]/track->fExpQ>3.5) continue;
1819     if (l<2 || l>3){
1820       deltan = (cl->GetNy()+cl->GetNz()-ny[l]-nz[l]);
1821     }
1822     else{
1823       deltan = (cl->GetNz()-nz[l]);
1824     }
1825     if (deltan>2.0) continue;  // extended - highly probable shared cluster
1826     //
1827     for (Int_t itrack=3;itrack>=0;itrack--){
1828       if (fgLayers[l].fClusterTracks[itrack][c]<0) continue;
1829       if (fgLayers[l].fClusterTracks[itrack][c]!=trackID){
1830        tracks[trackindex]  = fgLayers[l].fClusterTracks[itrack][c];
1831        trackindex++;
1832       }
1833     }
1834   }
1835   if (trackindex==0) return -1;
1836   if (trackindex==1){
1837     sharedtrack = tracks[0];
1838   }else{
1839     if (trackindex==2) sharedtrack =TMath::Min(tracks[0],tracks[1]);
1840     else{
1841       //
1842       Int_t track[24], cluster[24];
1843       for (Int_t i=0;i<trackindex;i++){ track[i]=-1; cluster[i]=0;}
1844       Int_t index =0;
1845       //
1846       for (Int_t i=0;i<trackindex;i++){
1847         if (tracks[i]<0) continue;
1848         track[index] = tracks[i];
1849         cluster[index]++;
1850         for (Int_t j=i+1;j<trackindex;j++){
1851           if (tracks[j]<0) continue;
1852           if (tracks[j]==tracks[i]){
1853             cluster[index]++;
1854             tracks[j]=-1;
1855           }
1856         }
1857         index++;
1858       }
1859       Int_t max=0;
1860       for (Int_t i=0;i<index;i++){
1861         if (cluster[index]>max) {
1862           sharedtrack=track[index];
1863           max=cluster[index];
1864         }
1865       }
1866     }
1867   }
1868
1869   if (sharedtrack>=100000) return -1;
1870   //
1871   // make list of overlaps
1872   shared =0;
1873   for (Int_t icluster=0;icluster<6;icluster++){
1874     if (clusterlist[icluster]<0) continue;
1875     Int_t index = clusterlist[icluster];
1876     Int_t l=(index & 0xf0000000) >> 28;
1877     Int_t c=(index & 0x0fffffff) >> 00;
1878     if (c>fgLayers[l].fN) continue;
1879     AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(index);
1880     if (l==0 || l==1){
1881       if (cl->GetNy()>2) continue;
1882       if (cl->GetNz()>2) continue;
1883     }
1884     if (l==4 || l==5){
1885       if (cl->GetNy()>3) continue;
1886       if (cl->GetNz()>3) continue;
1887     }
1888     //
1889     for (Int_t itrack=3;itrack>=0;itrack--){
1890       if (fgLayers[l].fClusterTracks[itrack][c]<0) continue;
1891       if (fgLayers[l].fClusterTracks[itrack][c]==sharedtrack){
1892         overlist[l]=index;
1893         shared++;
1894       }
1895     }
1896   }
1897   return sharedtrack;
1898 }
1899
1900
1901 AliITStrackMI *  AliITStrackerMI::GetBest2Tracks(Int_t trackID1, Int_t trackID2, Float_t th0, Float_t th1){
1902   //
1903   // try to find track hypothesys without conflicts
1904   // with minimal chi2;
1905   TClonesArray *arr1 = (TClonesArray*)fTrackHypothesys.At(trackID1);
1906   Int_t entries1 = arr1->GetEntriesFast();
1907   TClonesArray *arr2 = (TClonesArray*)fTrackHypothesys.At(trackID2);
1908   if (!arr2) return (AliITStrackMI*) arr1->UncheckedAt(0);
1909   Int_t entries2 = arr2->GetEntriesFast();
1910   if (entries2<=0) return (AliITStrackMI*) arr1->UncheckedAt(0);
1911   //
1912   AliITStrackMI * track10=(AliITStrackMI*) arr1->UncheckedAt(0);
1913   AliITStrackMI * track20=(AliITStrackMI*) arr2->UncheckedAt(0);
1914   if (TMath::Abs(1./track10->Get1Pt())>0.5+TMath::Abs(1/track20->Get1Pt())) return track10;
1915
1916   for (Int_t itrack=0;itrack<entries1;itrack++){
1917     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
1918     UnRegisterClusterTracks(track,trackID1);
1919   }
1920   //
1921   for (Int_t itrack=0;itrack<entries2;itrack++){
1922     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
1923     UnRegisterClusterTracks(track,trackID2);
1924   }
1925   Int_t index1=0;
1926   Int_t index2=0;
1927   Float_t maxconflicts=6;
1928   Double_t maxchi2 =1000.;
1929   //
1930   // get weight of hypothesys - using best hypothesy
1931   Double_t w1,w2;
1932
1933   Int_t list1[6],list2[6];
1934   AliITSclusterV2 *clist1[6], *clist2[6] ;
1935   RegisterClusterTracks(track10,trackID1);
1936   RegisterClusterTracks(track20,trackID2);
1937   Float_t conflict1 = GetNumberOfSharedClusters(track10,trackID1,list1,clist1);
1938   Float_t conflict2 = GetNumberOfSharedClusters(track20,trackID2,list2,clist2);
1939   UnRegisterClusterTracks(track10,trackID1);
1940   UnRegisterClusterTracks(track20,trackID2);
1941   //
1942   // normalized chi2
1943   Float_t chi21 =0,chi22=0,ncl1=0,ncl2=0;
1944   Float_t nerry[6],nerrz[6];
1945   Float_t *erry1=GetErrY(trackID1),*errz1 = GetErrZ(trackID1);
1946   Float_t *erry2=GetErrY(trackID2),*errz2 = GetErrZ(trackID2);
1947   for (Int_t i=0;i<6;i++){
1948      if ( (erry1[i]>0) && (erry2[i]>0)) {
1949        nerry[i] = TMath::Min(erry1[i],erry2[i]);
1950        nerrz[i] = TMath::Min(errz1[i],errz2[i]);
1951      }else{
1952        nerry[i] = TMath::Max(erry1[i],erry2[i]);
1953        nerrz[i] = TMath::Max(errz1[i],errz2[i]);
1954      }
1955      if (TMath::Abs(track10->fDy[i])>0.000000000000001){
1956        chi21 += track10->fDy[i]*track10->fDy[i]/(nerry[i]*nerry[i]);
1957        chi21 += track10->fDz[i]*track10->fDz[i]/(nerrz[i]*nerrz[i]);
1958        ncl1++;
1959      }
1960      if (TMath::Abs(track20->fDy[i])>0.000000000000001){
1961        chi22 += track20->fDy[i]*track20->fDy[i]/(nerry[i]*nerry[i]);
1962        chi22 += track20->fDz[i]*track20->fDz[i]/(nerrz[i]*nerrz[i]);
1963        ncl2++;
1964      }
1965   }
1966   chi21/=ncl1;
1967   chi22/=ncl2;
1968   //
1969   //
1970   Float_t d1 = TMath::Sqrt(track10->fD[0]*track10->fD[0]+track10->fD[1]*track10->fD[1])+0.1;
1971   Float_t d2 = TMath::Sqrt(track20->fD[0]*track20->fD[0]+track20->fD[1]*track20->fD[1])+0.1;
1972   Float_t s1 = TMath::Sqrt(track10->GetSigmaY2()*track10->GetSigmaZ2());
1973   Float_t s2 = TMath::Sqrt(track20->GetSigmaY2()*track20->GetSigmaZ2());
1974   //
1975   w1 = (d2/(d1+d2)+ 2*s2/(s1+s2)+
1976         +s2/(s1+s2)*0.5*(chi22+2.)/(chi21+chi22+4.)
1977         +1.*TMath::Abs(1./track10->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
1978         );
1979   w2 = (d1/(d1+d2)+ 2*s1/(s1+s2)+
1980         s1/(s1+s2)*0.5*(chi21+2.)/(chi21+chi22+4.)
1981         +1.*TMath::Abs(1./track20->Get1Pt())/(TMath::Abs(1./track10->Get1Pt())+TMath::Abs(1./track20->Get1Pt()))
1982         );
1983
1984   Double_t sumw = w1+w2;
1985   w1/=sumw;
1986   w2/=sumw;
1987   if (w1<w2*0.5) {
1988     w1 = (d2+0.5)/(d1+d2+1.);
1989     w2 = (d1+0.5)/(d1+d2+1.);
1990   }
1991   //  Float_t maxmax       = w1*track10->fChi2MIP[0]+w2*track20->fChi2MIP[0]+w1*conflict1+w2*conflict2+1.;
1992   //Float_t maxconflicts0 = w1*conflict1+w2*conflict2;
1993   //
1994   // get pair of "best" hypothesys
1995   //
1996   Float_t * ny1 = GetNy(trackID1), * nz1 = GetNz(trackID1);
1997   Float_t * ny2 = GetNy(trackID2), * nz2 = GetNz(trackID2);
1998
1999   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
2000     AliITStrackMI * track1=(AliITStrackMI*) arr1->UncheckedAt(itrack1);
2001     //if (track1->fFakeRatio>0) continue;
2002     RegisterClusterTracks(track1,trackID1);
2003     for (Int_t itrack2=0;itrack2<entries2;itrack2++){
2004       AliITStrackMI * track2=(AliITStrackMI*) arr2->UncheckedAt(itrack2);
2005
2006       //      Float_t current = w1*track1->fChi2MIP[0]+w2*track2->fChi2MIP[0];
2007       //if (track2->fFakeRatio>0) continue;
2008       Float_t nskipped=0;
2009       RegisterClusterTracks(track2,trackID2);
2010       Int_t list1[6],list2[6];
2011       AliITSclusterV2 *clist1[6], *clist2[6] ;
2012       Float_t cconflict1 = GetNumberOfSharedClusters(track1,trackID1,list1,clist1);
2013       Float_t cconflict2 = GetNumberOfSharedClusters(track2,trackID2,list2,clist2);
2014       UnRegisterClusterTracks(track2,trackID2);
2015       //
2016       if (track1->fConstrain) nskipped+=w1*track1->fNSkipped;
2017       if (track2->fConstrain) nskipped+=w2*track2->fNSkipped;
2018       if (nskipped>0.5) continue;
2019       //
2020       //if ( w1*conflict1+w2*conflict2>maxconflicts0) continue;
2021       if (conflict1+1<cconflict1) continue;
2022       if (conflict2+1<cconflict2) continue;
2023       Float_t conflict=0;
2024       Float_t sumchi2=0;
2025       Float_t sum=0;
2026       for (Int_t i=0;i<6;i++){
2027         //
2028         Float_t c1 =0.; // conflict coeficients
2029         Float_t c2 =0.;
2030         if (clist1[i]&&clist2[i]){
2031           Float_t deltan = 0;
2032           if (i<2 || i>3){
2033             deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-TMath::Max(ny1[i],ny2[i])-TMath::Max(nz1[i],nz2[i]));
2034           }
2035           else{
2036             deltan = (clist1[i]->GetNz()-TMath::Max(nz1[i],nz2[i]));
2037           }
2038           c1  = 2./TMath::Max(3.+deltan,2.);
2039           c2  = 2./TMath::Max(3.+deltan,2.);
2040         }
2041         else{
2042           if (clist1[i]){
2043             Float_t deltan = 0;
2044             if (i<2 || i>3){
2045               deltan = (clist1[i]->GetNy()+clist1[i]->GetNz()-ny1[i]-nz1[i]);
2046             }
2047             else{
2048               deltan = (clist1[i]->GetNz()-nz1[i]);
2049             }
2050             c1  = 2./TMath::Max(3.+deltan,2.);
2051             c2  = 0;
2052           }
2053
2054           if (clist2[i]){
2055             Float_t deltan = 0;
2056             if (i<2 || i>3){
2057               deltan = (clist2[i]->GetNy()+clist2[i]->GetNz()-ny2[i]-nz2[i]);
2058             }
2059             else{
2060               deltan = (clist2[i]->GetNz()-nz2[i]);
2061             }
2062             c2  = 2./TMath::Max(3.+deltan,2.);
2063             c1  = 0;
2064           }
2065         }
2066         //
2067         Double_t chi21=0,chi22=0;
2068         if (TMath::Abs(track1->fDy[i])>0.) {
2069           chi21 = (track1->fDy[i]/track1->fSigmaY[i])*(track1->fDy[i]/track1->fSigmaY[i])+
2070             (track1->fDz[i]/track1->fSigmaZ[i])*(track1->fDz[i]/track1->fSigmaZ[i]);
2071           //chi21 = (track1->fDy[i]*track1->fDy[i])/(nerry[i]*nerry[i])+
2072           //  (track1->fDz[i]*track1->fDz[i])/(nerrz[i]*nerrz[i]);
2073         }else{
2074           if (TMath::Abs(track1->fSigmaY[i]>0.)) c1=1;
2075         }
2076         //
2077         if (TMath::Abs(track2->fDy[i])>0.) {
2078           chi22 = (track2->fDy[i]/track2->fSigmaY[i])*(track2->fDy[i]/track2->fSigmaY[i])+
2079             (track2->fDz[i]/track2->fSigmaZ[i])*(track2->fDz[i]/track2->fSigmaZ[i]);
2080           //chi22 = (track2->fDy[i]*track2->fDy[i])/(nerry[i]*nerry[i])+
2081           //  (track2->fDz[i]*track2->fDz[i])/(nerrz[i]*nerrz[i]);
2082         }
2083         else{
2084           if (TMath::Abs(track2->fSigmaY[i]>0.)) c2=1;
2085         }
2086         sumchi2+=w1*(1.+c1)*(1+c1)*(chi21+c1)+w2*(1.+c2)*(1+c2)*(chi22+c2);
2087         if (chi21>0) sum+=w1;
2088         if (chi22>0) sum+=w2;
2089         conflict+=(c1+c2);
2090       }
2091       Double_t norm = sum-w1*track1->fNSkipped-w2*track2->fNSkipped;
2092       if (norm<0) norm =1/(w1*track1->fNSkipped+w2*track2->fNSkipped);
2093       Double_t normchi2 = 2*conflict+sumchi2/sum;
2094       if ( normchi2 <maxchi2 ){
2095         index1 = itrack1;
2096         index2 = itrack2;
2097         maxconflicts = conflict;
2098         maxchi2 = normchi2;
2099       }
2100     }
2101     UnRegisterClusterTracks(track1,trackID1);
2102   }
2103   //
2104   //  if (maxconflicts<4 && maxchi2<th0){
2105   if (maxchi2<th0*2.){
2106     Float_t orig = track10->fFakeRatio*track10->GetNumberOfClusters();
2107     AliITStrackMI* track1=(AliITStrackMI*) arr1->UncheckedAt(index1);
2108     track1->fChi2MIP[5] = maxconflicts;
2109     track1->fChi2MIP[6] = maxchi2;
2110     track1->fChi2MIP[7] = 0.01+orig-(track1->fFakeRatio*track1->GetNumberOfClusters());
2111     //    track1->UpdateESDtrack(AliESDtrack::kITSin);
2112     track1->fChi2MIP[8] = index1;
2113     fBestTrackIndex[trackID1] =index1;
2114     UpdateESDtrack(track1, AliESDtrack::kITSin);
2115   }
2116   else if (track10->fChi2MIP[0]<th1){
2117     track10->fChi2MIP[5] = maxconflicts;
2118     track10->fChi2MIP[6] = maxchi2;
2119     //    track10->UpdateESDtrack(AliESDtrack::kITSin);
2120     UpdateESDtrack(track10,AliESDtrack::kITSin);
2121   }
2122
2123   for (Int_t itrack=0;itrack<entries1;itrack++){
2124     AliITStrackMI * track=(AliITStrackMI*) arr1->UncheckedAt(itrack);
2125     UnRegisterClusterTracks(track,trackID1);
2126   }
2127   //
2128   for (Int_t itrack=0;itrack<entries2;itrack++){
2129     AliITStrackMI * track=(AliITStrackMI*) arr2->UncheckedAt(itrack);
2130     UnRegisterClusterTracks(track,trackID2);
2131   }
2132
2133   if (track10->fConstrain&&track10->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track10->fChi2MIP[1]<kMaxChi2PerCluster[1]
2134       &&track10->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track10->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2135     //  if (track10->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track10->fChi2MIP[1]<kMaxChi2PerCluster[1]
2136   //    &&track10->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track10->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2137     RegisterClusterTracks(track10,trackID1);
2138   }
2139   if (track20->fConstrain&&track20->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track20->fChi2MIP[1]<kMaxChi2PerCluster[1]
2140       &&track20->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track20->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2141     //if (track20->fChi2MIP[0]<kMaxChi2PerCluster[0]&&track20->fChi2MIP[1]<kMaxChi2PerCluster[1]
2142     //  &&track20->fChi2MIP[2]<kMaxChi2PerCluster[2]&&track20->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2143     RegisterClusterTracks(track20,trackID2);
2144   }
2145   return track10;
2146
2147 }
2148
2149 void AliITStrackerMI::UseClusters(const AliKalmanTrack *t, Int_t from) const {
2150   //--------------------------------------------------------------------
2151   // This function marks clusters assigned to the track
2152   //--------------------------------------------------------------------
2153   AliTracker::UseClusters(t,from);
2154
2155   AliITSclusterV2 *c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(0));
2156   //if (c->GetQ()>2) c->Use();
2157   if (c->GetSigmaZ2()>0.1) c->Use();
2158   c=(AliITSclusterV2 *)GetCluster(t->GetClusterIndex(1));
2159   //if (c->GetQ()>2) c->Use();
2160   if (c->GetSigmaZ2()>0.1) c->Use();
2161
2162 }
2163
2164
2165 void AliITStrackerMI::AddTrackHypothesys(AliITStrackMI * track, Int_t esdindex)
2166 {
2167   //------------------------------------------------------------------
2168   // add track to the list of hypothesys
2169   //------------------------------------------------------------------
2170
2171   if (esdindex>=fTrackHypothesys.GetEntriesFast()) fTrackHypothesys.Expand(esdindex*2+10);
2172   //
2173   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2174   if (!array) {
2175     array = new TObjArray(10);
2176     fTrackHypothesys.AddAt(array,esdindex);
2177   }
2178   array->AddLast(track);
2179 }
2180
2181 void AliITStrackerMI::SortTrackHypothesys(Int_t esdindex, Int_t maxcut, Int_t mode)
2182 {
2183   //-------------------------------------------------------------------
2184   // compress array of track hypothesys
2185   // keep only maxsize best hypothesys
2186   //-------------------------------------------------------------------
2187   if (esdindex>fTrackHypothesys.GetEntriesFast()) return;
2188   if (! (fTrackHypothesys.At(esdindex)) ) return;
2189   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2190   Int_t entries = array->GetEntriesFast();
2191   //
2192   //- find preliminary besttrack as a reference
2193   Float_t minchi2=10000;
2194   Int_t maxn=0;
2195   AliITStrackMI * besttrack=0;
2196   for (Int_t itrack=0;itrack<array->GetEntriesFast();itrack++){
2197     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2198     if (!track) continue;
2199     Float_t chi2 = NormalizedChi2(track,0);
2200     //
2201     Int_t tpcLabel=track->fESDtrack->GetTPCLabel();
2202     track->SetLabel(tpcLabel);
2203     CookdEdx(track);
2204     track->fFakeRatio=1.;
2205     CookLabel(track,0.); //For comparison only
2206     //
2207     //if (chi2<kMaxChi2PerCluster[0]&&track->fFakeRatio==0){
2208     if (chi2<kMaxChi2PerCluster[0]){
2209       if (track->GetNumberOfClusters()<maxn) continue;
2210       maxn = track->GetNumberOfClusters();
2211       if (chi2<minchi2){
2212         minchi2=chi2;
2213         besttrack=track;
2214       }
2215     }
2216     else{
2217       if (track->fConstrain || track->fN>5){  //keep best short tracks - without vertex constrain
2218         delete array->RemoveAt(itrack);
2219       }
2220     }
2221   }
2222   if (!besttrack) return;
2223   //
2224   //
2225   //take errors of best track as a reference
2226   Float_t *erry = GetErrY(esdindex), *errz = GetErrZ(esdindex);
2227   Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
2228   for (Int_t i=0;i<6;i++) {
2229     if (besttrack->fClIndex[i]>0){
2230       erry[i] = besttrack->fSigmaY[i]; erry[i+6] = besttrack->fSigmaY[i+6];
2231       errz[i] = besttrack->fSigmaZ[i]; errz[i+6] = besttrack->fSigmaZ[i+6];
2232       ny[i]   = besttrack->fNy[i];
2233       nz[i]   = besttrack->fNz[i];
2234     }
2235   }
2236   //
2237   // calculate normalized chi2
2238   //
2239   Float_t * chi2        = new Float_t[entries];
2240   Int_t * index         = new Int_t[entries];
2241   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2242   for (Int_t itrack=0;itrack<entries;itrack++){
2243     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2244     if (track){
2245       track->fChi2MIP[0] = GetNormalizedChi2(track, mode);
2246       if (track->fChi2MIP[0]<kMaxChi2PerCluster[0])
2247         chi2[itrack] = track->fChi2MIP[0];
2248       else{
2249         if (track->fConstrain || track->fN>5){  //keep best short tracks - without vertex constrain
2250           delete array->RemoveAt(itrack);
2251         }
2252       }
2253     }
2254   }
2255   //
2256   TMath::Sort(entries,chi2,index,kFALSE);
2257   besttrack = (AliITStrackMI*)array->At(index[0]);
2258   if (besttrack&&besttrack->fChi2MIP[0]<kMaxChi2PerCluster[0]){
2259     for (Int_t i=0;i<6;i++){
2260       if (besttrack->fClIndex[i]>0){
2261         erry[i] = besttrack->fSigmaY[i]; erry[i+6] = besttrack->fSigmaY[i+6];
2262         errz[i] = besttrack->fSigmaZ[i]; erry[i+6] = besttrack->fSigmaY[i+6];
2263         ny[i]   = besttrack->fNy[i];
2264         nz[i]   = besttrack->fNz[i];
2265       }
2266     }
2267   }
2268   //
2269   // calculate one more time with updated normalized errors
2270   for (Int_t i=0;i<entries;i++) chi2[i] =10000;
2271   for (Int_t itrack=0;itrack<entries;itrack++){
2272     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2273     if (track){
2274       track->fChi2MIP[0] = GetNormalizedChi2(track,mode);
2275       if (track->fChi2MIP[0]<kMaxChi2PerCluster[0])
2276         chi2[itrack] = track->fChi2MIP[0]-0*(track->GetNumberOfClusters()+track->fNDeadZone);
2277       else
2278         {
2279           if (track->fConstrain || track->fN>5){  //keep best short tracks - without vertex constrain
2280             delete array->RemoveAt(itrack);
2281           }
2282         }
2283     }
2284   }
2285   entries = array->GetEntriesFast();
2286   //
2287   //
2288   if (entries>0){
2289     TObjArray * newarray = new TObjArray();
2290     TMath::Sort(entries,chi2,index,kFALSE);
2291     besttrack = (AliITStrackMI*)array->At(index[0]);
2292     if (besttrack){
2293       //
2294       for (Int_t i=0;i<6;i++){
2295         if (besttrack->fNz[i]>0&&besttrack->fNy[i]>0){
2296           erry[i] = besttrack->fSigmaY[i]; erry[i+6] = besttrack->fSigmaY[i+6];
2297           errz[i] = besttrack->fSigmaZ[i]; errz[i+6] = besttrack->fSigmaZ[i+6];
2298           ny[i]   = besttrack->fNy[i];
2299           nz[i]   = besttrack->fNz[i];
2300         }
2301       }
2302       besttrack->fChi2MIP[0] = GetNormalizedChi2(besttrack,mode);
2303       Float_t minchi2 = TMath::Min(besttrack->fChi2MIP[0]+5.+besttrack->fNUsed, double(kMaxChi2PerCluster[0]));
2304       Float_t minn = besttrack->GetNumberOfClusters()-3;
2305       Int_t accepted=0;
2306       for (Int_t i=0;i<entries;i++){
2307         AliITStrackMI * track = (AliITStrackMI*)array->At(index[i]);
2308         if (!track) continue;
2309         if (accepted>maxcut) break;
2310         track->fChi2MIP[0] = GetNormalizedChi2(track,mode);
2311         if (track->fConstrain || track->fN>5){  //keep best short tracks - without vertex constrain
2312           if (track->GetNumberOfClusters()<6 && (track->fChi2MIP[0]+track->fNUsed>minchi2)){
2313             delete array->RemoveAt(index[i]);
2314             continue;
2315           }
2316         }
2317         Bool_t shortbest = !track->fConstrain && track->fN<6;
2318         if ((track->fChi2MIP[0]+track->fNUsed<minchi2 && track->GetNumberOfClusters()>=minn) ||shortbest){
2319           if (!shortbest) accepted++;
2320           //
2321           newarray->AddLast(array->RemoveAt(index[i]));
2322           for (Int_t i=0;i<6;i++){
2323             if (nz[i]==0){
2324               erry[i] = track->fSigmaY[i]; erry[i+6] = track->fSigmaY[i+6];
2325               errz[i] = track->fSigmaZ[i]; errz[i]   = track->fSigmaZ[i+6];
2326               ny[i]   = track->fNy[i];
2327               nz[i]   = track->fNz[i];
2328             }
2329           }
2330         }
2331         else{
2332           delete array->RemoveAt(index[i]);
2333         }
2334       }
2335       array->Delete();
2336       delete fTrackHypothesys.RemoveAt(esdindex);
2337       fTrackHypothesys.AddAt(newarray,esdindex);
2338     }
2339     else{
2340       array->Delete();
2341       delete fTrackHypothesys.RemoveAt(esdindex);
2342     }
2343   }
2344   delete [] chi2;
2345   delete [] index;
2346 }
2347
2348
2349
2350 AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI * original, Int_t checkmax)
2351 {
2352   //-------------------------------------------------------------
2353   // try to find best hypothesy
2354   // currently - minimal chi2 of track+backpropagated track+matching to the tpc track
2355   //-------------------------------------------------------------
2356   if (fTrackHypothesys.GetEntriesFast()<=esdindex) return 0;
2357   TObjArray * array = (TObjArray*) fTrackHypothesys.At(esdindex);
2358   if (!array) return 0;
2359   Int_t entries = array->GetEntriesFast();
2360   if (!entries) return 0;
2361   Float_t minchi2 = 100000;
2362   AliITStrackMI * besttrack=0;
2363   //
2364   AliITStrackMI * backtrack    = new AliITStrackMI(*original);
2365   AliITStrackMI * forwardtrack = new AliITStrackMI(*original);
2366   Double_t xyzv[]={GetX(),GetY(),GetZ()};
2367   Double_t ersv[]={GetSigmaX()/3.,GetSigmaY()/3.,GetSigmaZ()/3.};
2368   //
2369   for (Int_t i=0;i<entries;i++){
2370     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
2371     if (!track) continue;
2372     Float_t sigmarfi,sigmaz;
2373     GetDCASigma(track,sigmarfi,sigmaz);
2374     track->fDnorm[0] = sigmarfi;
2375     track->fDnorm[1] = sigmaz;
2376     //
2377     track->fChi2MIP[1] = 1000000;
2378     track->fChi2MIP[2] = 1000000;
2379     track->fChi2MIP[3] = 1000000;
2380     //
2381     // backtrack
2382     backtrack = new(backtrack) AliITStrackMI(*track);
2383     if (track->fConstrain){
2384       if (!backtrack->PropagateTo(3.,0.0028,65.19)) continue;
2385       if (!backtrack->Improve(0,xyzv,ersv))         continue;
2386       if (!backtrack->PropagateTo(2.,0.0028,0))     continue;
2387       if (!backtrack->Improve(0,xyzv,ersv))         continue;
2388       if (!backtrack->PropagateTo(1.,0.0028,0))     continue;
2389       if (!backtrack->Improve(0,xyzv,ersv))         continue;
2390       if (!backtrack->PropagateToVertex())          continue;
2391       backtrack->ResetCovariance();
2392       if (!backtrack->Improve(0,xyzv,ersv))         continue;
2393     }else{
2394       backtrack->ResetCovariance();
2395     }
2396     backtrack->ResetClusters();
2397
2398     Double_t x = original->GetX();
2399     if (!RefitAt(x,backtrack,track)) continue;
2400     //
2401     track->fChi2MIP[1] = NormalizedChi2(backtrack,0);
2402     //for (Int_t i=2;i<6;i++){track->fDy[i]+=backtrack->fDy[i]; track->fDz[i]+=backtrack->fDz[i];}
2403     if (track->fChi2MIP[1]>kMaxChi2PerCluster[1]*6.)  continue;
2404     track->fChi22 = GetMatchingChi2(backtrack,original);
2405
2406     if ((track->fConstrain) && track->fChi22>90.)  continue;
2407     if ((!track->fConstrain) && track->fChi22>30.)  continue;
2408     if ( track->fChi22/track->GetNumberOfClusters()>11.)  continue;
2409
2410
2411     if  (!(track->fConstrain)&&track->fChi2MIP[1]>kMaxChi2PerCluster[1])  continue;
2412     Bool_t isOK=kTRUE;
2413     if(!isOK) continue;
2414     //
2415     //forward track - without constraint
2416     forwardtrack = new(forwardtrack) AliITStrackMI(*original);
2417     forwardtrack->ResetClusters();
2418     x = track->GetX();
2419     RefitAt(x,forwardtrack,track);
2420     track->fChi2MIP[2] = NormalizedChi2(forwardtrack,0);
2421     if  (track->fChi2MIP[2]>kMaxChi2PerCluster[2]*6.0)  continue;
2422     if  (!(track->fConstrain)&&track->fChi2MIP[2]>kMaxChi2PerCluster[2])  continue;
2423
2424     track->fD[0] = forwardtrack->GetD(GetX(),GetY());
2425     track->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
2426     forwardtrack->fD[0] = track->fD[0];
2427     forwardtrack->fD[1] = track->fD[1];
2428     {
2429       Int_t list[6];
2430       AliITSclusterV2* clist[6];
2431       track->fChi2MIP[4] = GetNumberOfSharedClusters(track,esdindex,list,clist);
2432       if ( (!track->fConstrain) && track->fChi2MIP[4]>1.0) continue;
2433     }
2434
2435     track->fChi2MIP[3] = GetInterpolatedChi2(forwardtrack,backtrack);
2436     if  ( (track->fChi2MIP[3]>6.*kMaxChi2PerCluster[3])) continue;
2437     if  ( (!track->fConstrain) && (track->fChi2MIP[3]>2*kMaxChi2PerCluster[3])) {
2438       track->fChi2MIP[3]=1000;
2439       continue;
2440     }
2441     Double_t chi2 = track->fChi2MIP[0]+track->fNUsed;
2442     //
2443     for (Int_t ichi=0;ichi<5;ichi++){
2444       forwardtrack->fChi2MIP[ichi] = track->fChi2MIP[ichi];
2445     }
2446     if (chi2 < minchi2){
2447       //besttrack = new AliITStrackMI(*forwardtrack);
2448       besttrack = track;
2449       besttrack->SetLabel(track->GetLabel());
2450       besttrack->fFakeRatio = track->fFakeRatio;
2451       minchi2   = chi2;
2452       original->fD[0] = forwardtrack->GetD(GetX(),GetY());
2453       original->fD[1] = forwardtrack->GetZat(GetX())-GetZ();
2454     }
2455   }
2456   delete backtrack;
2457   delete forwardtrack;
2458   Int_t accepted=0;
2459   for (Int_t i=0;i<entries;i++){
2460     AliITStrackMI * track = (AliITStrackMI*)array->At(i);
2461     if (!track) continue;
2462
2463     if (accepted>checkmax || track->fChi2MIP[3]>kMaxChi2PerCluster[3]*6. ||
2464         (track->GetNumberOfClusters()<besttrack->GetNumberOfClusters()-1.)||
2465         track->fChi2MIP[0]>besttrack->fChi2MIP[0]+2.*besttrack->fNUsed+3.){
2466       if (track->fConstrain || track->fN>5){  //keep best short tracks - without vertex constrain
2467         delete array->RemoveAt(i);
2468         continue;
2469       }
2470     }
2471     else{
2472       accepted++;
2473     }
2474   }
2475   //
2476   array->Compress();
2477   SortTrackHypothesys(esdindex,checkmax,1);
2478   array = (TObjArray*) fTrackHypothesys.At(esdindex);
2479   besttrack = (AliITStrackMI*)array->At(0);
2480   if (!besttrack)  return 0;
2481   besttrack->fChi2MIP[8]=0;
2482   fBestTrackIndex[esdindex]=0;
2483   entries = array->GetEntriesFast();
2484   AliITStrackMI *longtrack =0;
2485   minchi2 =1000;
2486   Float_t minn=besttrack->GetNumberOfClusters()+besttrack->fNDeadZone;
2487   for (Int_t itrack=entries-1;itrack>0;itrack--){
2488     AliITStrackMI * track = (AliITStrackMI*)array->At(itrack);
2489     if (!track->fConstrain) continue;
2490     if (track->GetNumberOfClusters()+track->fNDeadZone<minn) continue;
2491     if (track->fChi2MIP[0]-besttrack->fChi2MIP[0]>0.0) continue;
2492     if (track->fChi2MIP[0]>4.) continue;
2493     minn = track->GetNumberOfClusters()+track->fNDeadZone;
2494     longtrack =track;
2495   }
2496   //if (longtrack) besttrack=longtrack;
2497
2498   Int_t list[6];
2499   AliITSclusterV2 * clist[6];
2500   Float_t shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
2501   if (besttrack->fConstrain&&besttrack->fChi2MIP[0]<kMaxChi2PerCluster[0]&&besttrack->fChi2MIP[1]<kMaxChi2PerCluster[1]
2502       &&besttrack->fChi2MIP[2]<kMaxChi2PerCluster[2]&&besttrack->fChi2MIP[3]<kMaxChi2PerCluster[3]){
2503     RegisterClusterTracks(besttrack,esdindex);
2504   }
2505   //
2506   //
2507   if (shared>0.0){
2508     Int_t nshared;
2509     Int_t overlist[6];
2510     Int_t sharedtrack = GetOverlapTrack(besttrack, esdindex, nshared, list, overlist);
2511     if (sharedtrack>=0){
2512       //
2513       besttrack = GetBest2Tracks(esdindex,sharedtrack,10,5.5);
2514       if (besttrack){
2515         shared = GetNumberOfSharedClusters(besttrack,esdindex,list,clist);
2516       }
2517       else return 0;
2518     }
2519   }
2520
2521   if (shared>2.5) return 0;
2522   if (shared>1.0) return besttrack;
2523   //
2524   // Don't sign clusters if not gold track
2525   //
2526   if (!besttrack->IsGoldPrimary()) return besttrack;
2527   if (besttrack->fESDtrack->GetKinkIndex(0)!=0) return besttrack;   //track belong to kink
2528   //
2529   if (fConstraint[fPass]){
2530     //
2531     // sign clusters
2532     //
2533     Float_t *ny = GetNy(esdindex), *nz = GetNz(esdindex);
2534     for (Int_t i=0;i<6;i++){
2535       Int_t index = besttrack->fClIndex[i];
2536       if (index<=0) continue;
2537       Int_t ilayer =  (index & 0xf0000000) >> 28;
2538       if (besttrack->fSigmaY[ilayer]<0.00000000001) continue;
2539       AliITSclusterV2 *c = (AliITSclusterV2*)GetCluster(index);
2540       if (!c) continue;
2541       if (ilayer>3&&c->GetNy()+c->GetNz()>6) continue;
2542       if ( (c->GetNy()+c->GetNz() )> ny[i]+nz[i]+0.7) continue; //shared track
2543       if (  c->GetNz()> nz[i]+0.7) continue; //shared track
2544       if ( ilayer>2&& besttrack->fNormQ[ilayer]/besttrack->fExpQ>1.5) continue;
2545       //if (  c->GetNy()> ny[i]+0.7) continue; //shared track
2546
2547       Bool_t cansign = kTRUE;
2548       for (Int_t itrack=0;itrack<entries; itrack++){
2549         AliITStrackMI * track = (AliITStrackMI*)array->At(i);
2550         if (!track) continue;
2551         if (track->fChi2MIP[0]>besttrack->fChi2MIP[0]+2.*shared+1.) break;
2552         if ( (track->fClIndex[ilayer]>0) && (track->fClIndex[ilayer]!=besttrack->fClIndex[ilayer])){
2553           cansign = kFALSE;
2554           break;
2555         }
2556       }
2557       if (cansign){
2558         if (TMath::Abs(besttrack->fDy[ilayer]/besttrack->fSigmaY[ilayer])>3.) continue;
2559         if (TMath::Abs(besttrack->fDz[ilayer]/besttrack->fSigmaZ[ilayer])>3.) continue;
2560         if (!c->IsUsed()) c->Use();
2561       }
2562     }
2563   }
2564   return besttrack;
2565 }
2566
2567
2568
2569 void  AliITStrackerMI::GetBestHypothesysMIP(TObjArray &itsTracks)
2570 {
2571   //
2572   // get "best" hypothesys
2573   //
2574
2575   Int_t nentries = itsTracks.GetEntriesFast();
2576   for (Int_t i=0;i<nentries;i++){
2577     AliITStrackMI* track = (AliITStrackMI*)itsTracks.At(i);
2578     if (!track) continue;
2579     TObjArray * array = (TObjArray*) fTrackHypothesys.At(i);
2580     if (!array) continue;
2581     if (array->GetEntriesFast()<=0) continue;
2582     //
2583     AliITStrackMI* longtrack=0;
2584     Float_t minn=0;
2585     Float_t maxchi2=1000;
2586     for (Int_t j=0;j<array->GetEntriesFast();j++){
2587       AliITStrackMI* track = (AliITStrackMI*)array->At(j);
2588       if (!track) continue;
2589       if (track->fGoldV0) {
2590         longtrack = track;   //gold V0 track taken
2591         break;
2592       }
2593       if (track->GetNumberOfClusters()+track->fNDeadZone<minn) continue;
2594       Float_t chi2 = track->fChi2MIP[0];
2595       if (fAfterV0){
2596         if (!track->fGoldV0&&track->fConstrain==kFALSE) chi2+=5;
2597       }
2598       if (track->GetNumberOfClusters()+track->fNDeadZone>minn) maxchi2 = track->fChi2MIP[0];
2599       //
2600       if (chi2 > maxchi2) continue;
2601       minn= track->GetNumberOfClusters()+track->fNDeadZone;
2602       maxchi2 = chi2;
2603       longtrack=track;
2604     }
2605     //
2606     //
2607     //
2608     AliITStrackMI * besttrack = (AliITStrackMI*)array->At(0);
2609     if (!longtrack) {longtrack = besttrack;}
2610     else besttrack= longtrack;
2611     //
2612     if (besttrack){
2613       Int_t list[6];
2614       AliITSclusterV2 * clist[6];
2615       Float_t shared = GetNumberOfSharedClusters(longtrack,i,list,clist);
2616       //
2617       track->fNUsed = shared;
2618       track->fNSkipped = besttrack->fNSkipped;
2619       track->fChi2MIP[0] = besttrack->fChi2MIP[0];
2620       if (shared>0){
2621         Int_t nshared;
2622         Int_t overlist[6];
2623         //
2624         Int_t sharedtrack = GetOverlapTrack(longtrack, i, nshared, list, overlist);
2625         //if (sharedtrack==-1) sharedtrack=0;
2626         if (sharedtrack>=0){
2627           besttrack = GetBest2Tracks(i,sharedtrack,10,5.5);
2628         }
2629       }
2630       if (besttrack&&fAfterV0){
2631         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
2632       }
2633       if (besttrack&&fConstraint[fPass])
2634         UpdateESDtrack(besttrack,AliESDtrack::kITSin);
2635       //if (besttrack&&besttrack->fConstrain)
2636       //        UpdateESDtrack(besttrack,AliESDtrack::kITSin);
2637       if (besttrack->fChi2MIP[0]+besttrack->fNUsed>1.5){
2638         if ( (TMath::Abs(besttrack->fD[0])>0.1) && fConstraint[fPass]) {
2639           track->fReconstructed= kFALSE;
2640         }
2641         if ( (TMath::Abs(besttrack->fD[1])>0.1) && fConstraint[fPass]){
2642           track->fReconstructed= kFALSE;
2643         }
2644       }
2645
2646     }
2647   }
2648 }
2649
2650
2651 void AliITStrackerMI::CookLabel(AliITStrackMI *track,Float_t wrong) const {
2652   //--------------------------------------------------------------------
2653   //This function "cooks" a track label. If label<0, this track is fake.
2654   //--------------------------------------------------------------------
2655   Int_t tpcLabel=-1;
2656
2657   if ( track->fESDtrack)   tpcLabel =  TMath::Abs(track->fESDtrack->GetTPCLabel());
2658
2659    track->fChi2MIP[9]=0;
2660    Int_t nwrong=0;
2661    for (Int_t i=0;i<track->GetNumberOfClusters();i++){
2662      Int_t cindex = track->GetClusterIndex(i);
2663      Int_t l=(cindex & 0xf0000000) >> 28;
2664      AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(cindex);
2665      Int_t isWrong=1;
2666      for (Int_t ind=0;ind<3;ind++){
2667        if (tpcLabel>0)
2668          if (cl->GetLabel(ind)==tpcLabel) isWrong=0;
2669      }
2670      track->fChi2MIP[9]+=isWrong*(2<<l);
2671      nwrong+=isWrong;
2672    }
2673    track->fFakeRatio = double(nwrong)/double(track->GetNumberOfClusters());
2674    if (tpcLabel>0){
2675      if (track->fFakeRatio>wrong) track->fLab = -tpcLabel;
2676      else
2677        track->fLab = tpcLabel;
2678    }
2679
2680 }
2681
2682
2683
2684 void AliITStrackerMI::CookdEdx(AliITStrackMI* track)
2685 {
2686   //
2687   //
2688   //  Int_t list[6];
2689   //AliITSclusterV2 * clist[6];
2690   //  Int_t shared = GetNumberOfSharedClusters(track,index,list,clist);
2691   Float_t dedx[4];
2692   Int_t accepted=0;
2693   track->fChi2MIP[9]=0;
2694   for (Int_t i=0;i<track->GetNumberOfClusters();i++){
2695     Int_t cindex = track->GetClusterIndex(i);
2696     Int_t l=(cindex & 0xf0000000) >> 28;
2697     AliITSclusterV2 *cl = (AliITSclusterV2*)GetCluster(cindex);
2698     Int_t lab = TMath::Abs(track->fESDtrack->GetTPCLabel());
2699     Int_t isWrong=1;
2700     for (Int_t ind=0;ind<3;ind++){
2701       if (cl->GetLabel(ind)==lab) isWrong=0;
2702     }
2703     track->fChi2MIP[9]+=isWrong*(2<<l);
2704     if (l<2) continue;
2705     //if (l>3 && (cl->GetNy()>4) || (cl->GetNz()>4)) continue;  //shared track
2706     //if (l>3&& !(cl->GetType()==1||cl->GetType()==10)) continue;
2707     //if (l<4&& !(cl->GetType()==1)) continue;
2708     dedx[accepted]= track->fdEdxSample[i];
2709     //dedx[accepted]= track->fNormQ[l];
2710     accepted++;
2711   }
2712   if (accepted<1) {
2713     track->SetdEdx(0);
2714     return;
2715   }
2716   Int_t indexes[4];
2717   TMath::Sort(accepted,dedx,indexes,kFALSE);
2718   Double_t low=0.;
2719   Double_t up=0.51;
2720   Double_t nl=low*accepted, nu =up*accepted;
2721   Float_t sumamp = 0;
2722   Float_t sumweight =0;
2723   for (Int_t i=0; i<accepted; i++) {
2724     Float_t weight =1;
2725     if (i<nl+0.1)     weight = TMath::Max(1.-(nl-i),0.);
2726     if (i>nu-1)     weight = TMath::Max(nu-i,0.);
2727     sumamp+= dedx[indexes[i]]*weight;
2728     sumweight+=weight;
2729   }
2730   track->SetdEdx(sumamp/sumweight);
2731 }
2732
2733
2734 void  AliITStrackerMI::MakeCoeficients(Int_t ntracks){
2735   //
2736   //
2737   if (fCoeficients) delete []fCoeficients;
2738   fCoeficients = new Float_t[ntracks*48];
2739   for (Int_t i=0;i<ntracks*48;i++) fCoeficients[i]=-1.;
2740 }
2741
2742
2743 Double_t AliITStrackerMI::GetPredictedChi2MI(AliITStrackMI* track, const AliITSclusterV2 *cluster,Int_t layer)
2744 {
2745   //
2746   //
2747   //
2748   Float_t erry,errz;
2749   Float_t theta = track->GetTgl();
2750   Float_t phi   = track->GetSnp();
2751   phi = TMath::Sqrt(phi*phi/(1.-phi*phi));
2752   GetError(layer,cluster,theta,phi,track->fExpQ,erry,errz);
2753   Double_t chi2 = track->GetPredictedChi2MI(cluster->GetY(),cluster->GetZ(),erry,errz);
2754   Float_t ny,nz;
2755   GetNTeor(layer,cluster, theta,phi,ny,nz);
2756   Double_t delta = cluster->GetNy()+cluster->GetNz()-nz-ny;
2757   if (delta>1){
2758     chi2+=0.5*TMath::Min(delta/2,2.);
2759     chi2+=2.*cluster->GetDeltaProbability();
2760   }
2761   //
2762   track->fNy[layer] =ny;
2763   track->fNz[layer] =nz;
2764   track->fSigmaY[layer] = erry;
2765   track->fSigmaZ[layer] = errz;
2766   //track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt(1+theta*theta+phi*phi);
2767   track->fNormQ[layer] = cluster->GetQ()/TMath::Sqrt((1.+ track->fP3*track->fP3)/(1.- track->fP2*track->fP2));
2768   return chi2;
2769
2770 }
2771
2772 Int_t    AliITStrackerMI::UpdateMI(AliITStrackMI* track, const AliITSclusterV2* cl,Double_t chi2,Int_t index) const
2773 {
2774   //
2775   //
2776   //
2777   Int_t layer = (index & 0xf0000000) >> 28;
2778   track->fClIndex[layer] = index;
2779   if ( (layer>1) &&track->fNormQ[layer]/track->fExpQ<0.5 ) {
2780     chi2+= (0.5-track->fNormQ[layer]/track->fExpQ)*10.;
2781     track->fdEdxMismatch+=(0.5-track->fNormQ[layer]/track->fExpQ)*10.;
2782   }
2783   return track->UpdateMI(cl->GetY(),cl->GetZ(),track->fSigmaY[layer],track->fSigmaZ[layer],chi2,index);
2784 }
2785
2786 void AliITStrackerMI::GetNTeor(Int_t layer, const AliITSclusterV2* /*cl*/, Float_t theta, Float_t phi, Float_t &ny, Float_t &nz)
2787 {
2788   //
2789   //get "mean shape"
2790   //
2791   if (layer==0){
2792     ny = 1.+TMath::Abs(phi)*3.2;
2793     nz = 1.+TMath::Abs(theta)*0.34;
2794     return;
2795   }
2796   if (layer==1){
2797     ny = 1.+TMath::Abs(phi)*3.2;
2798     nz = 1.+TMath::Abs(theta)*0.28;
2799     return;
2800   }
2801
2802   if (layer>3){
2803     ny = 2.02+TMath::Abs(phi)*1.95;
2804     nz = 2.02+TMath::Abs(phi)*2.35;
2805     return;
2806   }
2807   ny  = 6.6-2.7*TMath::Abs(phi);
2808   nz  = 2.8-3.11*TMath::Abs(phi)+0.45*TMath::Abs(theta);
2809 }
2810
2811
2812
2813 Int_t AliITStrackerMI::GetError(Int_t layer, const AliITSclusterV2*cl, Float_t theta, Float_t phi,Float_t expQ, Float_t &erry, Float_t &errz)
2814 {
2815   //calculate cluster position error
2816   //
2817   Float_t nz,ny;
2818   GetNTeor(layer, cl,theta,phi,ny,nz);
2819   erry   = TMath::Sqrt(cl->GetSigmaY2());
2820   errz   = TMath::Sqrt(cl->GetSigmaZ2());
2821   //
2822   // PIXELS
2823   if (layer<2){
2824
2825     if (TMath::Abs(ny-cl->GetNy())>0.6)  {
2826       if (ny<cl->GetNy()){
2827         erry*=0.4+TMath::Abs(ny-cl->GetNy());
2828         errz*=0.4+TMath::Abs(ny-cl->GetNy());
2829       }else{
2830         erry*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
2831         errz*=0.7+0.5*TMath::Abs(ny-cl->GetNy());
2832       }
2833     }
2834     if (TMath::Abs(nz-cl->GetNz())>1.)  {
2835       erry*=TMath::Abs(nz-cl->GetNz());
2836       errz*=TMath::Abs(nz-cl->GetNz());
2837     }
2838     erry*=0.85;
2839     errz*=0.85;
2840     erry= TMath::Min(erry,float(0.005));
2841     errz= TMath::Min(errz,float(0.03));
2842     return 10;
2843   }
2844
2845 //STRIPS
2846   if (layer>3){
2847     //factor 1.8 appears in new simulation
2848     //
2849     Float_t scale=1.8;
2850     if (cl->GetNy()==100||cl->GetNz()==100){
2851       erry = 0.004*scale;
2852       errz = 0.2*scale;
2853       return 100;
2854     }
2855     if (cl->GetNy()+cl->GetNz()>12){
2856       erry = 0.06*scale;
2857       errz = 0.57*scale;
2858       return 100;
2859     }
2860     Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
2861     Float_t chargematch = TMath::Max(double(normq/expQ),2.);
2862     //
2863     if (cl->GetType()==1 || cl->GetType()==10 ){
2864       if (chargematch<1.0 || (cl->GetNy()+cl->GetNz()<nz+ny+0.5)){
2865         errz = 0.043*scale;
2866         erry = 0.00094*scale;
2867         return 101;
2868       }
2869       if (cl->GetNy()+cl->GetNz()<nz+ny+1.2){
2870         errz = 0.06*scale;
2871         erry =0.0013*scale;
2872         return 102;
2873       }
2874       erry = 0.0027*scale;
2875       errz = TMath::Min(0.028*(chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.15)*scale;
2876       return 103;
2877     }
2878     if (cl->GetType()==2 || cl->GetType()==11 ){
2879       erry = TMath::Min(0.0010*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.05)*scale;
2880       errz = TMath::Min(0.025*(1+chargematch+cl->GetNy()+cl->GetNz()-nz+ny),0.5)*scale;
2881       return 104;
2882     }
2883
2884     if (cl->GetType()>100 ){
2885       if ((chargematch+cl->GetNy()+cl->GetNz()-nz-ny<1.5)){
2886         errz = 0.05*scale;
2887         erry = 0.00096*scale;
2888         return 105;
2889       }
2890       if (cl->GetNy()+cl->GetNz()-nz-ny<1){
2891         errz = 0.10*scale;
2892         erry = 0.0025*scale;
2893         return 106;
2894       }
2895
2896       errz = TMath::Min(0.05*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.4)*scale;
2897       erry = TMath::Min(0.003*(chargematch+cl->GetNy()+cl->GetNz()-nz-ny),0.05)*scale;
2898       return 107;
2899     }
2900     Float_t diff = cl->GetNy()+cl->GetNz()-ny-nz;
2901     if (diff<1) diff=1;
2902     if (diff>4) diff=4;
2903
2904     if (cl->GetType()==5||cl->GetType()==6||cl->GetType()==7||cl->GetType()==8){
2905       errz = 0.14*diff;
2906       erry = 0.003*diff;
2907       return 108;
2908     }
2909     erry = 0.04*diff;
2910     errz = 0.06*diff;
2911     return 109;
2912   }
2913   //DRIFTS
2914   Float_t normq = cl->GetQ()/(TMath::Sqrt(1+theta*theta+phi*phi));
2915   Float_t chargematch = normq/expQ;
2916   Float_t factorz=1;
2917   Int_t   cnz = cl->GetNz()%10;
2918   //charge match
2919   if (cl->GetType()==1){
2920     if (chargematch<1.25){
2921       erry =  0.0028*(1.+6./cl->GetQ());  // gold clusters
2922     }
2923     else{
2924       erry = 0.003*chargematch;
2925       if (cl->GetNz()==3) erry*=1.5;
2926     }
2927     if (chargematch<1.0){
2928       errz =  0.0011*(1.+6./cl->GetQ());
2929     }
2930     else{
2931       errz = 0.002*(1+2*(chargematch-1.));
2932     }
2933     if (cnz>nz+0.6) {
2934       erry*=(cnz-nz+0.5);
2935       errz*=1.4*(cnz-nz+0.5);
2936     }
2937   }
2938   if (cl->GetType()>1){
2939     if (chargematch<1){
2940       erry =  0.00385*(1.+6./cl->GetQ());  // gold clusters
2941       errz =  0.0016*(1.+6./cl->GetQ());
2942     }
2943     else{
2944       errz = 0.0014*(1+3*(chargematch-1.));
2945       erry = 0.003*(1+3*(chargematch-1.));
2946     }
2947     if (cnz>nz+0.6) {
2948       erry*=(cnz-nz+0.5);
2949       errz*=1.4*(cnz-nz+0.5);
2950     }
2951   }
2952
2953   if (TMath::Abs(cl->GetY())>2.5){
2954     factorz*=1+2*(TMath::Abs(cl->GetY())-2.5);
2955   }
2956   if (TMath::Abs(cl->GetY())<1){
2957     factorz*=1.+0.5*TMath::Abs(TMath::Abs(cl->GetY())-1.);
2958   }
2959   factorz= TMath::Min(factorz,float(4.));
2960   errz*=factorz;
2961
2962   erry= TMath::Min(erry,float(0.05));
2963   errz= TMath::Min(errz,float(0.05));
2964   return 200;
2965 }
2966
2967
2968
2969 void   AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Float_t &sigmaz)
2970 {
2971   //
2972   //DCA sigmas parameterization
2973   //to be paramterized using external parameters in future
2974   //
2975   //
2976   sigmarfi = 0.004+1.4 *TMath::Abs(track->fP4)+332.*track->fP4*track->fP4;
2977   sigmaz   = 0.011+4.37*TMath::Abs(track->fP4);
2978
2979 }
2980
2981
2982 void AliITStrackerMI::SignDeltas( TObjArray *ClusterArray, Float_t vz)
2983 {
2984   //
2985   //
2986   Int_t entries = ClusterArray->GetEntriesFast();
2987   if (entries<4) return;
2988   AliITSclusterV2* cluster = (AliITSclusterV2*)ClusterArray->At(0);
2989   Int_t layer = cluster->GetLayer();
2990   if (layer>1) return;
2991   Int_t index[10000];
2992   Int_t ncandidates=0;
2993   Float_t r = (layer>0)? 7:4;
2994   //
2995   for (Int_t i=0;i<entries;i++){
2996     AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(i);
2997     Float_t nz = 1+TMath::Abs((cl0->GetZ()-vz)/r);
2998     if (cl0->GetNy()+cl0->GetNz()<=5+2*layer+nz) continue;
2999     index[ncandidates] = i;  //candidate to belong to delta electron track
3000     ncandidates++;
3001     if (cl0->GetNy()+cl0->GetNz()>9+2*layer+nz) {
3002       cl0->SetDeltaProbability(1);
3003     }
3004   }
3005   //
3006   //
3007   //
3008   for (Int_t i=0;i<ncandidates;i++){
3009     AliITSclusterV2* cl0 = (AliITSclusterV2*)ClusterArray->At(index[i]);
3010     if (cl0->GetDeltaProbability()>0.8) continue;
3011     //
3012     Int_t ncl = 0;
3013     Float_t y[100],z[100],sumy,sumz,sumy2, sumyz, sumw;
3014     sumy=sumz=sumy2=sumyz=sumw=0.0;
3015     for (Int_t j=0;j<ncandidates;j++){
3016       if (i==j) continue;
3017       AliITSclusterV2* cl1 = (AliITSclusterV2*)ClusterArray->At(index[j]);
3018       //
3019       Float_t dz = cl0->GetZ()-cl1->GetZ();
3020       Float_t dy = cl0->GetY()-cl1->GetY();
3021       if (TMath::Sqrt(dz*dz+dy*dy)<0.2){
3022         Float_t weight = cl1->GetNy()+cl1->GetNz()-2;
3023         y[ncl] = cl1->GetY();
3024         z[ncl] = cl1->GetZ();
3025         sumy+= y[ncl]*weight;
3026         sumz+= z[ncl]*weight;
3027         sumy2+=y[ncl]*y[ncl]*weight;
3028         sumyz+=y[ncl]*z[ncl]*weight;
3029         sumw+=weight;
3030         ncl++;
3031       }
3032     }
3033     if (ncl<4) continue;
3034     Float_t det = sumw*sumy2  - sumy*sumy;
3035     Float_t delta=1000;
3036     if (TMath::Abs(det)>0.01){
3037       Float_t z0  = (sumy2*sumz - sumy*sumyz)/det;
3038       Float_t k   = (sumyz*sumw - sumy*sumz)/det;
3039       delta = TMath::Abs(cl0->GetZ()-(z0+k*cl0->GetY()));
3040     }
3041     else{
3042       Float_t z0  = sumyz/sumy;
3043       delta = TMath::Abs(cl0->GetZ()-z0);
3044     }
3045     if ( delta<0.05) {
3046       cl0->SetDeltaProbability(1-20.*delta);
3047     }
3048   }
3049 }
3050
3051
3052 void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
3053 {
3054   //
3055   //
3056   track->UpdateESDtrack(flags);
3057   AliITStrackMI * oldtrack = (AliITStrackMI*)(track->fESDtrack->GetITStrack());
3058   if (oldtrack) delete oldtrack;
3059   track->fESDtrack->SetITStrack(new AliITStrackMI(*track));
3060   if (TMath::Abs(track->fDnorm[1])<0.000000001){
3061     printf("Problem\n");
3062   }
3063 }
3064
3065
3066
3067
3068
3069
3070
3071
3072 void  AliITStrackerMI::FindV0(AliESD *event)
3073 {
3074   //
3075   // fast V0 finder
3076   //
3077   //TTreeSRedirector cstream("itsv0.root");
3078   Int_t centries=0;
3079   AliHelix helixes[30000];
3080   TObjArray trackarray(30000);
3081   TObjArray trackarrayc(30000);
3082   Float_t dist[30000];
3083   Float_t normdist0[30000];
3084   Float_t normdist1[30000];
3085   Float_t normdist[30000];
3086   Float_t norm[30000];
3087   AliESDV0MI  *vertexarray    = new AliESDV0MI[100000];
3088   AliESDV0MI *pvertex     = &vertexarray[0];
3089   AliITStrackMI * dummy=0;
3090   //
3091   //
3092   Int_t entries = fTrackHypothesys.GetEntriesFast();
3093   for (Int_t i=0;i<entries;i++){
3094     TObjArray * array = (TObjArray*)fTrackHypothesys.At(i);
3095     if (!array) continue;
3096     // get best track without vertex constrain
3097     Int_t hentries = array->GetEntriesFast();
3098     //
3099     // best with vertex constrain
3100     AliITStrackMI * trackc = (AliITStrackMI*)array->At(0);
3101     if (trackc&&trackc->fConstrain&&trackc->fN==6&&trackc->fNormChi2[0]<2.) continue;
3102     trackc=0;
3103     for (Int_t ih=0;ih<hentries;ih++){
3104       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3105       if (!trackh->fConstrain) continue;
3106       if (trackh->fN<6) continue;
3107       trackc = trackh;
3108       if (!dummy) dummy = trackc;
3109       dummy->SetLabel(0);
3110       break;
3111     }
3112     //
3113     // best without vertex
3114     AliITStrackMI * track = 0;
3115     for (Int_t ih=0;ih<hentries;ih++){
3116       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
3117       if (trackh->fConstrain) continue;
3118       track = trackh;
3119       break;
3120     }
3121     if (trackc&&track){
3122       if (trackc->fChi2MIP[1]<2.) continue;
3123       if (trackc->fChi2MIP[0]<2. && trackc->fChi2MIP[1]<2.) continue;
3124       trackarrayc.AddAt(trackc,i);
3125       if (trackc->fN==6&&track->fN&&trackc->fNormChi2[0] < track->fNormChi2[0]-2) continue;
3126     }
3127     //
3128     //
3129     //
3130     if (track){
3131       dist[i] = TMath::Sqrt(track->fD[0]*track->fD[0]+track->fD[1]*track->fD[1]);
3132       norm[i] = track->fDnorm[0];
3133       normdist0[i] = TMath::Abs(track->fD[0]/track->fDnorm[0]);
3134       normdist1[i] = TMath::Abs(track->fD[1]/track->fDnorm[1]);
3135       normdist[i]  = TMath::Sqrt(normdist0[i]*normdist0[i]+normdist1[i]*normdist1[i]);
3136       if (track->IsGoldPrimary()) continue; //primary track
3137       if (track->fD[0]<0.02 && (track->fN+track->fNDeadZone>5.8)){
3138         if (normdist[i]<3.) continue;  // primary track  - cutoff 3 sigma
3139         if (normdist0[i]<2.) continue; //DCA normalized cut 2 sigma
3140       }
3141       trackarray.AddAt(track,i);
3142       new (&helixes[i]) AliHelix(*track);
3143     }
3144   }
3145   //
3146   //
3147   //  Int_t multifound=0;
3148   Int_t vertexall =0;
3149   AliESDV0MI tempvertex;
3150   Float_t primvertex[3]={GetX(),GetY(),GetZ()};
3151
3152
3153   for (Int_t itrack0=0;itrack0<entries;itrack0++){
3154     //
3155     AliITStrackMI * track0 = (AliITStrackMI*)trackarray.At(itrack0);
3156     if (!track0) continue;
3157     if (track0->fP4>0) continue;
3158     AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
3159     //
3160     TObjArray * array0     = (TObjArray*)fTrackHypothesys.At(itrack0);
3161     //
3162     Int_t vertexes =0;
3163     for (Int_t itrack1=0;itrack1<entries;itrack1++){
3164       AliITStrackMI * track1 = (AliITStrackMI*)trackarray.At(itrack1);
3165       if (!track1) continue;
3166       if (track1->fP4<0) continue;
3167       AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
3168       if (trackc0&&trackc1){
3169         if (TMath::Min(trackc0->fChi2MIP[1],trackc1->fChi2MIP[1])<2.) continue;
3170       }
3171       if (track1->fNDeadZone+track0->fNDeadZone>1.1) continue;
3172       TObjArray * array1     = (TObjArray*)fTrackHypothesys.At(itrack1);
3173       //
3174       //if (normdist0[itrack0]+normdist0[itrack1]<3) continue;
3175       //if (normdist[itrack0]+normdist[itrack1]<4) continue;
3176       //
3177       //
3178       AliHelix *h1 = &helixes[itrack0];
3179       AliHelix *h2 = &helixes[itrack1];
3180       Double_t rmin =0;
3181       Double_t distance = TestV0(h1,h2,pvertex,rmin);
3182       //
3183       if (distance>0.4) continue;
3184       if (pvertex->GetRr()<0.3) continue;
3185       if (pvertex->GetRr()>20.) continue;
3186       pvertex->SetM(*track0);
3187       pvertex->SetP(*track1);
3188       pvertex->Update(primvertex);
3189       if (pvertex->GetRr()<0.3) continue;
3190       if (pvertex->GetRr()>20.) continue;
3191       if (track1->fNDeadZone+track0->fNDeadZone>0.5 &&distance>0.12) continue;
3192
3193       //
3194
3195       if ( TMath::Abs((TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel())))<2
3196            ||(centries<5000&&(pvertex->GetPointAngle()>0.95))){
3197         //cstream<<"Iter0"<<track0<<track1<<pvertex<<normdist[itrack0]<<normdist[itrack1]<<"\n";
3198         centries++;
3199       }
3200       //
3201       //
3202       if (pvertex->GetPointAngle()<0.85) continue;
3203       //      if (normdist[itrack0]+normdist[itrack1]<6&&pvertex->GetPointAngle()<0.99) continue;
3204       //
3205       //
3206       pvertex->SetLab(0,track0->GetLabel());
3207       pvertex->SetLab(1,track1->GetLabel());
3208       pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
3209       pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
3210       // calculate chi2s
3211       //
3212       pvertex->SetChi2After(0);
3213       pvertex->SetChi2Before(0);
3214       pvertex->SetNBefore(0);
3215       pvertex->SetNAfter(0);
3216
3217       for (Int_t i=0;i<6;i++){
3218         Float_t radius = fgLayers[i].GetR();
3219         if (pvertex->GetRr()>radius+0.5){
3220           pvertex->SetNBefore(pvertex->GetNBefore()+2.);
3221           //
3222           if (track0->fClIndex[i]<=0) {
3223             pvertex->SetChi2Before(pvertex->GetChi2Before()+9);
3224           }else{
3225             Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
3226               track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
3227             pvertex->SetChi2Before(pvertex->GetChi2Before()+chi2);
3228           }
3229
3230           if (track1->fClIndex[i]<=0) {
3231             pvertex->SetChi2Before(pvertex->GetChi2Before()+9);
3232
3233           }else{
3234             Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
3235               track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
3236             //        pvertex->fChi2Before+=chi2;
3237             pvertex->SetChi2Before(pvertex->GetChi2Before()+chi2);
3238           }
3239         }
3240
3241         if (pvertex->GetRr()<radius-0.5){
3242           pvertex->SetNAfter(pvertex->GetNAfter()+2.);
3243           //
3244           if (track0->fClIndex[i]<=0) {
3245             pvertex->SetChi2After(pvertex->GetChi2After()+9);
3246           }else{
3247             Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
3248               track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
3249             pvertex->SetChi2After(pvertex->GetChi2After()+chi2);
3250           }
3251
3252           if (track1->fClIndex[i]<=0) {
3253             pvertex->SetChi2After(pvertex->GetChi2After()+9.);
3254           }else{
3255             Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
3256               track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
3257             pvertex->SetChi2After(pvertex->GetChi2After()+chi2);
3258           }
3259         }
3260       }
3261       if (pvertex->GetNBefore()>2){
3262         if (pvertex->GetChi2Before()/pvertex->GetNBefore()<4.) continue; //have clusters before vetex
3263       }
3264       Int_t ibest0=0,ibest1=0;
3265       AliITStrackMI * ntrack0 = track0;
3266       AliITStrackMI * ntrack1 = track1;
3267       //
3268       //
3269       //PH      Float_t oldistance = pvertex->GetDist2();
3270       Bool_t improve  = FindBestPair(itrack0,itrack1,pvertex,ibest0,ibest1);   // try to improve vertex
3271       if (pvertex->GetPointAngle()<0.5) continue;
3272       distance = pvertex->GetDist2();
3273       if (improve){
3274         ntrack0 = (AliITStrackMI*)array0->At(ibest0);
3275         ntrack1 = (AliITStrackMI*)array1->At(ibest1);
3276       }
3277       Bool_t accept0 = kFALSE;      // accept ==>  because of pointing angle
3278       if (pvertex->GetPointAngle()>0.999){
3279         if (pvertex->GetRr()<3.5 && (ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone)<11.5) continue;
3280         if (pvertex->GetRr()>3.5&& pvertex->GetDistNorm()<12) accept0 = kTRUE;
3281         if (pvertex->GetRr()>1 && pvertex->GetDist2()<0.1 && pvertex->GetDistNorm()<12) accept0 = kTRUE;
3282         if (pvertex->GetPointAngle()>0.9995&&pvertex->GetRr()>5.) accept0 = kTRUE;
3283       }
3284       Bool_t reject1= kFALSE;      // reject ==> bad kinematic
3285       //
3286       reject1 |=  TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>1.02 ||
3287                    TMath::Abs(ntrack0->fN-ntrack1->fN)>1.02;  // cut1
3288       reject1 |=  ntrack0->fNUsed+ntrack1->fNUsed>1.01;                                 // cut2
3289       reject1 |=  pvertex->GetDistNorm()>12;                                                // cut3
3290       reject1 |=  pvertex->GetDist2()>0.1 && improve;                                       // cut4
3291       reject1 |=  (TMath::Abs(ntrack0->fD[0])+TMath::Abs(ntrack1->fD[0]))/pvertex->GetDist2()<5; //cut5
3292       reject1 |=  TMath::Abs(ntrack0->fD[0]/pvertex->GetDist2())<2 || TMath::Abs(ntrack1->fD[0]/pvertex->GetDist2())<2;  //cut 6
3293       //
3294       // small radii cuts
3295       Bool_t reject2 = kFALSE;
3296       if (pvertex->GetRr()<3.6){
3297         reject2 |=  (TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>0.01);  // cut7
3298         reject2 |=  ntrack0->fNUsed+ntrack1->fNUsed>0.01;                                  // cut8
3299         reject2 |=  ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone<11.5;  // cut9
3300         reject2 |=  (ntrack0->fN+ntrack1->fN<11.5)&&pvertex->GetRr()<2;                        // cut10
3301         reject2 |=  pvertex->GetDist2()>0.1;                                                   // cut11
3302       }
3303       //PH      AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
3304       //PH      AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
3305       //
3306       //
3307       //
3308       if ( TMath::Abs((TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel())))<2
3309            ||(centries<500000)){
3310         /*
3311         cstream<<"It1"<<"Tr0.="<<ntrack0<<"TR1.="<<ntrack1<<"V0.="<<pvertex<<"ND0.="<<normdist[itrack0]<<"ND1.="<<
3312           normdist[itrack1]<<"D.="<<distance<<"DistOld="<<oldistance<<"Imp.="<<improve<<
3313           "A0="<<accept0<<"R1="<<reject1<<"R2="<<reject2<<"Rmin.="<<rmin<<
3314           "TrC0.="<<htrackc0<<"TRC1.="<<htrackc1<<"\n";
3315         */
3316         centries++;
3317       }
3318
3319       if (!accept0 && (reject1 || reject2))  continue;
3320
3321 //       if (distance>0.5) continue;
3322 //       distance = pvertex->GetDist2();
3323 //       if (pvertex->GetRr()>25 || pvertex->GetRr()<0.2) continue;
3324 //       if (pvertex->GetRr()/pvertex->fDistSigma<1) continue;
3325 //       if (pvertex->GetDistNorm()>10) continue;
3326 //       if (pvertex->GetPointAngle()<0.85) continue;
3327 //       if ((normdist[itrack0]<3||normdist[itrack1]<3)){
3328 //      if (pvertex->GetPointAngle()<0.99||pvertex->GetDist2()>0.15) continue;
3329 //       }
3330 //       if (distance>0.05*(0.8+0.2*(0.5+pvertex->GetRr()))) continue;
3331 //       if (pvertex->GetRr()<0.3) continue;
3332 //       if (pvertex->GetRr()>27.) continue;
3333
3334
3335
3336       //
3337       if (distance<0.3 &&pvertex->GetPointAngle()>0.998){
3338         track0->fGoldV0 = kTRUE;
3339         track1->fGoldV0 = kTRUE;
3340       }
3341       vertexes++;
3342       vertexall++;
3343       if (vertexall>=100000) break;
3344       pvertex = &vertexarray[vertexall];
3345     }
3346   }
3347   //  printf("\n\n\nMultifound\t%d\n\n\n",multifound);
3348   //
3349   // sort vertices according quality
3350   Float_t quality[40000];
3351   Int_t   indexes[40000];
3352   Int_t   trackvertices[40000];
3353   for (Int_t i=0;i<entries;i++) trackvertices[i]=0;
3354   for (Int_t i=0;i<vertexall;i++) {
3355     Float_t norm     = 1.-0.999*TMath::Abs(vertexarray[i].GetPointAngle());
3356     Float_t fnormdist = 1./(1+vertexarray[i].GetRr());
3357     quality[i] = norm*fnormdist;
3358   }
3359   //
3360   TMath::Sort(vertexall,quality,indexes,kFALSE);
3361
3362   for (Int_t i=0;i<vertexall;i++){
3363     pvertex= &vertexarray[indexes[i]];
3364     Int_t index0 = vertexarray[indexes[i]].GetIndex(0);
3365     Int_t index1 = vertexarray[indexes[i]].GetIndex(1);
3366     vertexarray[indexes[i]].SetOrder(2,i);
3367     vertexarray[indexes[i]].SetOrder(1,trackvertices[index1]);
3368     vertexarray[indexes[i]].SetOrder(0,trackvertices[index0]);
3369     Int_t v0index = event->AddV0MI(&vertexarray[indexes[i]]);
3370     //
3371     if (trackvertices[index1]+trackvertices[index0]>5) continue;
3372     if (trackvertices[index0]>2) continue;
3373     if (trackvertices[index1]>2) continue;
3374
3375     if (trackvertices[index1]+trackvertices[index0]>0) {
3376       //      if (pvertex->GetPointAngle()<0.995)  continue;
3377     }
3378     trackvertices[index0]++;
3379     trackvertices[index1]++;
3380
3381     AliESDtrack * ptrack0 = event->GetTrack(vertexarray[indexes[i]].GetIndex(0));
3382     AliESDtrack * ptrack1 = event->GetTrack(vertexarray[indexes[i]].GetIndex(1));
3383     if (!ptrack0 || !ptrack1){
3384       printf("BBBBBBBUUUUUUUUUUGGGGGGGGGG\n");
3385     }
3386     Int_t v0index0[3]={ptrack0->GetV0Index(0),ptrack0->GetV0Index(1),ptrack0->GetV0Index(2)};
3387     Int_t v0index1[3]={ptrack1->GetV0Index(0),ptrack1->GetV0Index(1),ptrack1->GetV0Index(2)};
3388     for (Int_t i=0;i<3;i++){
3389       if (v0index0[i]<0) {
3390         v0index0[i]=v0index;
3391         ptrack0->SetV0Indexes(v0index0);
3392         break;
3393       }
3394     }
3395     for (Int_t i=0;i<3;i++){
3396       if (v0index1[i]<0) {
3397         v0index1[i]=v0index;
3398         ptrack1->SetV0Indexes(v0index1);
3399         break;
3400       }
3401     }
3402   }
3403   delete[] vertexarray;
3404 }
3405
3406
3407
3408 Bool_t  AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1, AliESDV0MI *vertex, Int_t &i0, Int_t &i1)
3409 {
3410   //
3411   // try to find best pair from the tree of track hyp.
3412   //
3413   TObjArray * array0 = (TObjArray*)fTrackHypothesys.At(esdtrack0);
3414   Int_t entries0 = array0->GetEntriesFast();
3415   TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(esdtrack1);
3416   Int_t entries1 = array1->GetEntriesFast();
3417   //  AliITStrackMI *orig0 = (AliITStrackMI*)fOriginal.At(esdtrack0);
3418   //AliITStrackMI *orig1 = (AliITStrackMI*)fOriginal.At(esdtrack1);
3419   Double_t criticalradius = vertex->GetRr();
3420   AliITStrackMI * track0= 0;
3421   AliITStrackMI * track1= 0;
3422   i0 = -1;
3423   i1 = -1;
3424   //
3425   //
3426   Float_t rfirst0[2000];  //radius position of  the first cluster - track0
3427   Float_t rfirst1[2000];  //                                      - track1
3428   Float_t maxlocalx0=0;      //local x  for first  track
3429   Float_t maxlocalx1=0;      //local x  for second track
3430   Float_t cs0=1, sn0=0;            //rotations
3431   Float_t cs1=1, sn1=0;            //rotations
3432
3433   //
3434   for (Int_t itrack0=0;itrack0<entries0;itrack0++){
3435     rfirst0[itrack0]=-1.;
3436     AliITStrackMI * htrack0 = (AliITStrackMI*)array0->At(itrack0);
3437     if (!htrack0) continue;
3438     if (htrack0->fConstrain) continue;
3439     if (i0<0){
3440       i0     = itrack0;
3441       track0 = htrack0;
3442     }
3443     Double_t cs = TMath::Cos(htrack0->fAlpha);
3444     Double_t sn = TMath::Sin(htrack0->fAlpha);
3445     Double_t x  = htrack0->fX*cs - htrack0->fP0*sn;
3446     Double_t y  = htrack0->fX*sn + htrack0->fP0*cs;
3447     Double_t radius = TMath::Sqrt(x*x+y*y);
3448     if (criticalradius<3  && radius>6&&htrack0->fNDeadZone<0.2) continue;  // all cluster required
3449     if (criticalradius>10 && radius<6) continue;  // causality
3450     Double_t localx =  TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn);
3451     if (localx>maxlocalx0) {
3452       maxlocalx0=localx;
3453       cs0 = cs; sn0=sn;
3454     }
3455     rfirst0[itrack0] = radius;
3456   }
3457   for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3458     rfirst1[itrack1]=-1.;
3459     AliITStrackMI * htrack1 = (AliITStrackMI*)array1->At(itrack1);
3460     if (!htrack1) continue;
3461     if (htrack1->fConstrain) continue;
3462     if (i1<0){
3463       i1     = itrack1;
3464       track1 = htrack1;
3465     }
3466     Double_t cs = TMath::Cos(htrack1->fAlpha);
3467     Double_t sn = TMath::Sin(htrack1->fAlpha);
3468     //
3469     //
3470     Double_t x  = htrack1->fX*cs - htrack1->fP0*sn;
3471     Double_t y  = htrack1->fX*sn + htrack1->fP0*cs;
3472     Double_t radius = TMath::Sqrt(x*x+y*y);
3473     if (criticalradius<3  && radius>6&&htrack1->fNDeadZone<0.2) continue; //all clusters required
3474     if (criticalradius>10 && radius<6) continue; //causality
3475     Double_t localx =  TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn);
3476     if (localx>maxlocalx1) {
3477       maxlocalx1=localx;
3478       cs1 = cs; sn1=sn;
3479     }
3480     rfirst1[itrack1] = radius;
3481   }
3482   //
3483   //
3484   //
3485   const Float_t radiuses[4]={4,6.5,15.03,24.};
3486   //
3487   //
3488   // find the best tracks after decay point
3489   Float_t bestquality =100000;
3490   Float_t bestradius=0;
3491   AliESDV0MI v0;
3492   Double_t vpos[3];
3493   Float_t v[3]={GetX(),GetY(),GetZ()};
3494   //
3495   //
3496   for (Int_t itrack0=0;itrack0<entries0;itrack0++){
3497     if (rfirst0[itrack0]<0) continue;
3498     AliITStrackMI * htrack0 = (AliITStrackMI*)array0->At(itrack0);
3499     if (!htrack0) continue;
3500     //
3501     for (Int_t itrack1=0;itrack1<entries1;itrack1++){
3502       if (rfirst1[itrack1]<0) continue;
3503       AliITStrackMI * htrack1 = (AliITStrackMI*)array1->At(itrack1);
3504       if (!htrack1) continue;
3505       if (TMath::Abs(rfirst0[itrack0]-rfirst1[itrack1])>1.) continue;
3506       if (htrack0->fClIndex[6-htrack0->fN]==htrack1->fClIndex[6-htrack1->fN]) continue; //sharing of last cluster not allowe
3507       //
3508       if (htrack0->fNUsed+htrack1->fNUsed>0) continue; //sharing of clusters not allowed
3509       //
3510       //
3511       v0.SetM(*htrack0);
3512       v0.SetP(*htrack1);
3513       //      if (v0.Update(v)==0) continue;
3514       v0.Update(v);
3515       if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) <v0.GetRr()-0.3) continue;
3516       //
3517       //
3518       if (v0.GetDist2()>0.3) continue;
3519       if (v0.GetRr()<radiuses[1] && ( TMath::Abs(htrack0->fN+htrack0->fNDeadZone-htrack1->fN-htrack1->fNDeadZone)>0.5)) continue;
3520       //
3521       //if (v0.GetRr()<3. && (htrack0->fN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<11.7) continue;
3522       //if (v0.GetRr()<6. && (htrack0->fN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<9.7) continue;
3523       if (v0.GetRr()<3. && (htrack0->fN+htrack1->fN)<11.7) continue;
3524       if (v0.GetRr()<6. && (htrack0->fN+htrack1->fN)<9.7) continue;
3525       Double_t localx0=v0.GetXr(0)*cs0+v0.GetXr(1)*sn0;
3526       Double_t localx1=v0.GetXr(0)*cs1+v0.GetXr(1)*sn1;
3527       Double_t maxlocalx = TMath::Max(localx0,localx1);
3528       if (maxlocalx<3.4  && (htrack0->fN+htrack1->fN)<11.7) continue;
3529       if (maxlocalx<6.1  && (htrack0->fN+htrack1->fN)<9.7) continue;
3530       //
3531       Float_t fnormdist = v0.GetDist2()/0.05;
3532       fnormdist +=(htrack0->fNormChi2[6-htrack0->fN]+htrack1->fNormChi2[6-htrack1->fN]);
3533       fnormdist +=3*(htrack0->fNUsed+htrack1->fNUsed);
3534       if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) <v0.GetRr()+0.2){
3535         fnormdist +=(v0.GetRr()+0.2-TMath::Min(rfirst0[itrack0],rfirst1[itrack1]))/0.1;
3536       }
3537       //
3538       Float_t quality = fnormdist;
3539       if (quality<bestquality && v0.GetDist2()<vertex->GetDist2()){
3540         i0=itrack0;
3541         i1=itrack1;
3542         track0 =htrack0;
3543         track1 =htrack1;
3544         bestquality = quality;
3545         bestradius  = v0.GetRr();
3546         vpos[0] = v0.GetXr(0);
3547         vpos[1] = v0.GetXr(1);
3548         vpos[2] = v0.GetXr(2);
3549       }
3550     }
3551   }
3552   //
3553   //
3554   if (!track0||!track1) return kFALSE;
3555   if (track0->fNUsed+track1->fNUsed>0) return kFALSE;  // sharing of clusters not allowed
3556   //
3557   //propagate to vertex
3558
3559   Double_t alpha  = TMath::ATan2(vpos[1],vpos[0]);
3560   //
3561   AliITStrackMI track0p = *track0;
3562   AliITStrackMI track1p = *track1;
3563   //
3564   //  RefitAt(bestradius+0.5,&track0p,track0);
3565   //RefitAt(bestradius+0.5,&track1p,track1);
3566
3567   if ( track0p.Propagate(alpha,bestradius+0.2)) {
3568     track0= &track0p;
3569   }
3570   if (track1p.Propagate(alpha,bestradius+0.2)){
3571     track1 = &track1p;
3572   }
3573   //
3574   v0.SetM(*track0);
3575   v0.SetP(*track1);
3576   v0.Update(v);
3577   if (v0.GetDist2()<vertex->GetDist2() && v0.GetRr()<20){
3578     vertex->SetM(*track0);
3579     vertex->SetP(*track1);
3580     vertex->Update(v);
3581     return kTRUE;
3582   }
3583   return kFALSE;
3584 }
3585
3586
3587
3588
3589 Double_t   AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliESDV0MI *vertex, Double_t &rmin)
3590 {
3591   //
3592   // test the helixes for the distnce calculate vertex characteristic
3593   //
3594   rmin =0;
3595   Float_t distance1,distance2;
3596   AliHelix & dhelix1 = *helix1;
3597   Double_t pp[3],xx[3];
3598   dhelix1.GetMomentum(0,pp,0);
3599   dhelix1.Evaluate(0,xx);
3600   AliHelix &mhelix = *helix2;
3601   //
3602   //find intersection linear
3603   //
3604   Double_t phase[2][2],radius[2];
3605   Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
3606   Double_t delta1=10000,delta2=10000;
3607
3608   if (points>0){
3609     rmin = radius[0];
3610     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
3611     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
3612     dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
3613   }
3614   if (points==2){
3615     if (radius[1]<rmin) rmin = radius[1];
3616     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
3617     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
3618     dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
3619   }
3620   rmin = TMath::Sqrt(rmin);
3621   distance1 = TMath::Min(delta1,delta2);
3622   vertex->SetDist1(TMath::Sqrt(distance1));
3623
3624   //
3625   //find intersection parabolic
3626   //
3627   points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
3628   delta1=10000,delta2=10000;
3629
3630   if (points>0){
3631     dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
3632   }
3633   if (points==2){
3634     dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
3635   }
3636
3637   distance2 = TMath::Min(delta1,delta2);
3638   vertex->SetDist2(TMath::Sqrt(distance2));
3639   if (distance2<100){
3640     if (delta1<delta2){
3641       //get V0 info
3642       dhelix1.Evaluate(phase[0][0],vertex->GetXrp());
3643       dhelix1.GetMomentum(phase[0][0],vertex->GetPPp());
3644       mhelix.GetMomentum(phase[0][1],vertex->GetPMp());
3645       dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],vertex->GetAnglep());
3646       vertex->SetRr(TMath::Sqrt(radius[0]));
3647     }
3648     else{
3649       dhelix1.Evaluate(phase[1][0],vertex->GetXrp());
3650       dhelix1.GetMomentum(phase[1][0],vertex->GetPPp());
3651       mhelix.GetMomentum(phase[1][1],vertex->GetPMp());
3652       dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],vertex->GetAnglep());
3653       vertex->SetRr(TMath::Sqrt(radius[1]));
3654     }
3655   }
3656   //
3657   //
3658   return  vertex->GetDist2();
3659 }
3660
3661
3662
3663
3664
3665
3666