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