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