ITS/AliITStrackerMI.cxx

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