ITS/AliITStrackerMI.cxx

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