ITS/AliITStrackerMI.cxx

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